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This is a final technical report for the period 1 April 1978 - 
30 April 1979 for the NASA Grant NSG-7329 "Application of Laser 
Ranging and VLBI Data to a Study of Plate Tectonic Driving Forces." 

As outlined in our original proposal, our broad goals have been 
to investigate the conditions under which changes in plate driving 
or resistive forces associated with plate boundary earthquakes will 
be measurable with laser ranging or VLBI and to pinpoint those 
aspects of plate forces that can be characterized by such measurements. 

The bulk of this technical report consists of detailed reports 
of three aspects of our work as of the end of the grant year: 

(1) analytic solutions for two-dimensional stress diffusion in 
a plate following earthquake faulting on a finite fault; (2) 
two-dimensional finite-element solutions for the global state of 
stress at the Earth's surface for possible plate driving forces; 
and (3) finite-element solutions for three-dimensional stress 
diffusion in a viscoelastic Earth following earthquake faulting. 
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Introduction 

The migration of large earthquakes along plate boundaries 
has been a significant discovery in global tectonics (Mogi, 
1968; Kelleher, 1970; Delsemme and Smith, 1979). Since the 
occurrence of plate boundary earthquakes is associated with 
the state of stress at the convergingplate boundaries , it 
is necessary to understand the response of the lithosphere 
to application or removal of horizontal boundary stresses 
(or displacements) in space and time. The new geodetic 
data offered by the VLBI network will give us the opportunity 
to observe these displacements in time. 

There have been several studies on . stress and strain 
diffusion from plate boundaries . The results depend largely on 
the type of material properties used to represent the 
asthenosphere . Bolt and Dean (1973) investigated a one 
dimensional model of an elastic plate underlain by a viscous 
layer. Th^ confirm that the viscous shearing stress of the 
asthenosphere on the base of the lithosphere plays an important 
role in stress and strain diffusion through the lithosphere. 
They also predict the possible occurrence ' of a new type of 
strain wave in the lithosphere. 
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Melosh (1976) extends the study to a nonlinear viscous 

as thenosphere , again without an elastic component. Unlike 

the previous two viscous models, we introduce a Maxwellian 

behavior to the as thenosphere for the one-dimensional model; 

thus, an initial elastic response is possible. The results 

are significantly different from the previous models and 

suggest that, with the thickness of the asthenosphere at 

20 

about 80 km and a viscosity of 10 poise, the fluctuations 
in velocity and stress are greater than 50% of the equilibrium 
level for distances within 1000 km from the boundary '-disturbance 
(final report, NASA grant no. NSG-7329, April 1978). 

One-dimensional models, however, can provide only limited 
insight into the plate response; the- model ignores the 
finite length of the fault rupture. Moreover, the migration 
of major earthquakes is essentially parallel to the plate 
boundaries. In this study we examine displacements and 
stress propagation in two directions, parallel and perpendicular 
to the plate boundaries. Using this simple model will allow 
a qualitative understanding of the effect introduced by a 
finite fault on the plate dynamics . 
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Two dimensional stress and strain propagation . 

The ea.rth's plates can be considered as shells covering 
the spherical surface of the earth. In order to simplify 
the problem, the earth's curvature is neglected and plane 
stress is assumed within the lithosphere. Figure 1 shows 
a finite two dimensional plate of thickness overlying a 
second layer, the asthenosphere , of thickness The 

boundary x=0 corresponds to the converging boundary (or trench) 
and the boundary x=L represents the ridge. The width of the 
plate in the y direction (along the arc) is 2C. Now consider 
a small element of the plate •with volume dV=dxdydz. 

The forces in the x and y directions are 



The forces along the lithosphere per unit are are 



For the driving mechanism of the plates , the above forces 
must be balanced by interplate friction and basal shear 
stresses. In most studies the basal shear stress is a 
primary factor in balancing the forces (Richter, 1977; Davies, 
1978), the interplate friction appear to be smaller (Davies, 
1978). Although it is still not known whether the basal 
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shear drives or drags the plates, we will assume that it 
drags the plate motion and ignore the transform friction. 

If the as thenosphere is a simple elastic media with 
shear modulus yU , the basal shear on the top of the 
asthenosphere per unit area is 




( 2 ) 
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where "U. and IT are displacements in x and y directions 
respectively, and ^ is the thickness of the. asthenosphere. 

Combini/g (1) and (2) . and writing the . stresses in terms' of 
displacements, we obtain the equations of motion: 
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We can further simplify the equations of motion into 
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where 'v/ 7 is the Laplacian operator in two dimensions, and 
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Solving for and X <j ) gives the displacement fields 

u (x,y) and v(x,y). 

A set of solutions that satisfies (5), (6), (7), and (8) 
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The coefficients A^, A 2 , are functions of m. 

Boundary conditions 

In order to obtain non-trivial solutions, the boundary 
conditions require careful attention. The problems of 
geometry restrict the applications to general cases with 
rectangular boundaries. Great earthquakes and their aftershocks 
only occury a small portion of the boundary. To simulate this 
boundary condition, we assume a unit displacement along the 
x=0 corresponding to the arc. The displacement is. applied 
from y=-d to y=d throughout the thickness H^. The remainder 
of this boundary remains fixed. l-le assume frictionless boundarie 
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at y=-C . The ridge (x=L) is also assumed to be stress free. 


Mathematically , these conditions become 
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Using these boundary conditions we obtain: 
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Solving for A 2 and is accomplished with the matrix equation: 
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The solution has the form of a Fourier series. In practice 
we truncate the series to obtain a solution. The truncation 
introduces a Gibb's phenomena which .is reduced through a 
Lanczo’s factor. With 120 terms , the solutions for the 
displacements and stresses at the boundary x=0 (the 'arc') 
are plotted in Figure 2. The. figure shows u displacement 
applied as the boundary condition on the converging 'arc'. 

The corresponding displacement v along the 'arc' is shown in 
Figure 3. The displacements increase sharply near the "earthquake”. 
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Examining the stress distribution along the boundary 
in Figure 4, we note that the stress normal to the arc ( Cy ) 
is compressional in the 'earthquake' region and changes to 
tension at the edge of the rupture zone. 

The distribution of stress along the arc x=0 

(figure 5) fits our interpretation of the displacements in 
figure 3. A sharp increase in shear stress at the edge 
of the ’earthquake 1 rupture is also expected and is shown in 
figure 6. In each case these are the elastic deformations 
expected for the simple plate model at the plate boundary. 

We now wish to examine how deeply this instantaneous 
elastic deformation penetrates into the plate. Figure 7 and 8 
are the dis tributicns of the displacement u and stress 
along the x axis, that is, perpendicular to the arc. Both 
diagrams suggest that the elastic deformations are significant 
within a distance equivalent to the length of the fault 
rupture. This qualitatively agrees with exact solutions for 
the elastic deformations. If we now transform our elastic 
solution with the correspondence principle of viscoelasticity, 
we can obtain the time-dependent solutions for a Maxwellian 
as thenosphere . Our approach is analogous to the one- 
dimensional problem outlined in the previous report (final 
report, NASA grant no. NSG-7329). 
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FIGURES 

Figure 1 

Two dimensional model for an elastic lithosphere 
of thickness and Young's modulus E overlying a 
viscoelastic as thenosphere of thickness , viscosity , 
and shear modulusyW . The lithosphere has a finite length L 
in the x direction perdendicular to the arc. The width is 
2C in the y direction along the 'arc'. A unit displacement 
simulates the 'earthquake' at the boundary x=0 between 
y=-d and y=d. Consequently^ the earthquake rupture zone has 
length 2d. 

Figure 2 

Displacement u (in x direction) at the boundary x=0 
corresponding to the 'arc'. The displacement is plotted 
along the y axis, distance along the 'arc'. For this model 
we assume that d = 800 km and C = 4000 km. The length 
of the plate is 5000 km . The displacement u along the 
rupture is 10 m, 


Figure 3 

Displacement v parallel to the arc along the boundary x=0 
In the 'rupture' zone, the displacement increases rapidly 
until the termination of the rupture is reached at y= T d. 

Beyond the rupture, the displacement v decreases. 
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Figure 4 • 

The stress 0^ normal to the arc as a function of 
distance along the arc (x=0) . As expected, compression occurs 
in the earthquake zone and tension occurs outside this area. 

If the converging boundary is initially under compression, 
the 'earthquake' rupture increases the stress level outside 
the earthquake zone . 

Figure 5 

The stress parallel to the arc as a function of 

distance along the boundary (x=0) . 

< 

Figure 6 

Shear stress as a function of distance along the arc 
(x=0). Near the edge of the rupture zone the shear stress 
reaches a maximum. 

Figure 7 

Displacement normal to the arc (u) as a function of x, 
the distance from the converging arc; 

Figure 8 . 

Stress G~ x normal to the 'arc' as a function of x, the 
distance from the converging arc. 
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INTRODUCTION 

The state of stress in the lithosphere is the product 
of the many forces on a variety of scales that 
act on the plates. Foremost in importance among such 
forces are those responsible for plate motions, often 
collectively termed the driving mechanism. For specific 
models of the plate tectonic driving mechanism, the intra- 
plate stress field may be calculated and compared against 
those indicators of stress in the earth most likely to 
reflect tectonic* forces. The comparison of predicted 
and observed stresses, a powerful test of driving force 
models, is the subject of this paper. 

In earlier works, we calculated intraplate 
stresses for various driving mechanism models, using a 
finite difference technique, a coarse grid for the earth, 
and a body of observations consisting primarily of midplate 
earthquake mechanisms [ Solomon et al ., 1975,1977b; Richard - 
son et al . , 1976]. In this paper we offer several substantial 
improvements over the earlier work, including the use of a 
more versatile finite element computational technique with 
better spatial resolution, particularly near plate boundaries; 
a more extensive catalog of intraplate stress observations; 
and a wider range of tested driving force models. 

We begin with a review of the important tectonic driving 
and resistive forces capable of stressing the plates, with 
due consideration also for sources of intraplate stress not 
related to the driving mechanism. The constraints set on 
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these forces by the relative and absolute velocities of 
the plates are briefly summarized. 

The various methods for inferring the orientations 
and/or magnitudes of the principal intralithospheric stresses 
are next reviewed, including earthquake mechanisms, in 
situ measurements, and observation of stress-sensitive 
geological features. Based on these techniques, the 
regional stress orientation field and its observational 
basis is synthesized for each plate. 

The logical leap to the notion that the directions of 

* 

principal stresses inferred from observations may be used . 

as quantitative tests of driving force models as posed in 

this paper is based on three premises:. (1) that a regionally 

consistent stress field exists within large fractions of the 

lithosphere; (2) that such a stress field is steady on a 

6 

time period of less than'vlO years; and (3) that the domi- 
nant forces which give rise to most of that stress field 
are plate tectonic in origin. 

The first premise, regional consistency of stress 
orientations, is the most readily established for several 
important intraplate areas, on the basis of many deter- 
minations by a variety of • independent techniques . The 
second premise, steadiness of the regional stress field 
over time periods long compared to the period of observation, 
depends on the response of the plate to plate-boundary earthquakes 
through the rheology of the asthenosphere and lithosphere 
[c.f., Anderson , 1975] and is not well established at this 
time. Determination, by laser ranging [Bender and Silverberg, 


1975] or by very long baseline interferometry [Coates 



et al . , 1975] that the motion is steady between the stable 
interiors of different plates would lend strong support 
to the concept of a steady regional stress field. 

The third premise, that an identifiable portion of 
the regional stress field arises from plate tectonic 
forces, is the most difficult to ascertain. We argue 
below that this premise is most likely to be correct in 
regions well removed from such other notable sources of 
stress as recent tectonic or thermal activity, recent 
topographic loading or unloading, and pronounced struc- 
tural heterogeneities. 

If these three premises hold, then the comparison 
of calculated to observed intraplate stresses offers a 
sensitive test of driving mechanism models. If one or 
more of the premises is unfounded, then the .calculations 
of stresses presented in this paper may still be viewed 
as a guide to that contribution to the stress field from 
steady plate driving forces. 

The stress calculations are based on a finite element 
technique, described below and also in Richardson [1978] . 

The lithosphere is modeled as a thin shell, and the plates 
are regarded as distinct elastic units which terminate at 
plate boundaries. Net driving forces at spreading centers 
and subduction zones are abstracted as line forces at plate 
boundaries. Global intraplate stresses determined for a 
large number of possible force models, and the comparison 
of predictions with observations, are used as a basis to 
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modify and reject models. 

The results of the modeling may be summarized as follows: 
pushing forces related to the topography of ridges are 
required of all models that match the inferred intraplate 
stress field. The net pulling force exerted on the plates 
by subducted lithosphere is at most a few times greater than 
other forces acting on the plates, even if net slab forces 
decrease with increasing subduction rate because of resis- 
tive ' forces that depend on velocity. Forces acting at 
continental convergence zones that resist further convergence 
are important for models of the intraplate stress field 
in Europe, Asia, and the Indian plate. Viscous drag at the 
ba.se of the plate is best modeled as a resistive force, 
with a drag coefficient that is non-zero beneath oceanic 
lithosphere, but which may be concentrated by a factor of 
up to ten beneath continental lithosphere. Models of the 
driving mechanism in which drag forces drive the plates 
are in poor agreement with the observed data for most 
regions. A model in which drag balances the net torque 
produced on each plate by boundary forces, including slab 
pull concentrated on the subducted plate, yields good agree- 
ment between predicted and observed stresses in the contin- 
ents but a poor fit in oceanic regions. 
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DRIVING FORCES AND STRESSES 


The specification of the driving mechanism for plate 
tectonics should ideally be a fully dynamic description, in 
which plate motions and plate geometries do not have to be 
assumed but follow from a quantitative theory for mantle 
convection that includes plate recycling as a natural element. 

Such a specification is unfortunately not yet within our 
grasp. As an alternative, perhaps adequate at this stage in 
our knowledge, we abstract the true driving, mechanism as simply 
parameterized forces acting at plate boundaries and at the 
base of the lithosphere. The plate boundary forces depend 
upon boundary type, and perhaps on other parameters such as 
relative velocity or local plate age, and the basal forces 
depend on the nature of the interaction between the lithosphere 
and the asthenosphere . . 

Ridge Forces 

Mid-ocean ridges can exert a net compressive stress on 
the adjacent spreading plates because of their elevation with 
respect to abyssal sea floor. The magnitude of this stress 
has been estimated at several hundred bars [ Hales , 1969; Frank , 
1972; McKenzie , 1972; Artyushkov , 1973], and likely dominates 
any resistance to spreading at the ridge axis [ Sleep and Biehler , 
1970; Lachenbruch , 1973] as expressed in the extensional tectonics 
[Sykes , 1967] . 

The driving force associated with spreading centers is 
distributed over oceanic lithosphere young enough so that the 
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plate density structure continues to change with age [e.g., 

Hager , 1978], but we approximate it for simplicity as a line 
force at the ridge. The force is normal to the topographic 
strike of the ridge (i.e., generally in the direction of 
spreading) , and is symmetric about the ridge axis. Since 
the height of most ridges far from hot spots is nearly 
uniform [ Sclater et al ., 1971], the ridge force is approximately 
independent of spreading rate. 

Subduction Zone Forces 

A major potential driving force. for plate motions is the 
negative buoyancy of subducted lithosphere. The density 
contrast due to the temperature difference between subducted 
slab and surrounding mantle is capable of exerting a tensile 
stress of up to several kilobars on the subducted surface 
plate [ McKenzie , 1969a; Turcotte and Schubert , 1971], an 
important component of which is the elevation of major phase 
changes within the slab [ Smith and Toksoz , 1972; Solomon and 
Paw U , 1975]. 

For the determination of intraplate stress, the effective 
force at subduction zones is the net difference between 
the negative buoyancy of the slab and the variety of resistive 
forces acting to oppose subduction. These resistive forces 
include friction along the island arc thrust zone, drag against 
the sides of the slab, and the excess pressure on the leading 
edge of the slab [ Isacks and Molnar , 1971; Smith and Toksoz , 
1972; Richter , 1977]. It is also possible that large normal- 
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faulting earthquakes may decouple sections of the slab from the 
surface plate [ Kanamori , 1971; Abe , 1972]. Evidence both from 
the speeds of subducted plates and from the focal mechanisms 
of intraslab earthquakes suggest that the negative slab 
buoyancy is approximately balanced by these resistive forces 

II 

for most subduction zones [ Smith and Toksoz , 1972; Forsyth 
and Uyeda , 1975; Richter , 1977]. The net pulling force 
at trenches resulting from this approximate balance is likely 
to depend on convergence rate, the angle of subduction, and 
the age of the subducted lithosphere [ Richter and McKenzie , 

1979]. 

The extent to which the subducted slab exerts a 'puli' 
on the overthrust plate is an important but poorly resolved 
issue. That such a force exists is strongly suggested by 
the trenchward component of absolute motion for nearly all 
overthrust sides of subduction zones, particularly if 
inter-arc spreading is considered [ Chase , 1978b] . This force 
may be associated with a, convective cell generated by shear on 
the asthenosphere by the" descending slab [ Richter , 1973] . 

The relative magnitudes of net pull on the subducted plate 
and any pull on the overriding plate are not well constrained. 

The topography of the outer rise seaward of many trenches 
has been explained in terms of plate flexure. Hanks [1971] 
argued on the basis of a thin elastic plate model that both 
large (kilobars ) bending stresses and a superposed regional 
compressive stress of comparable magnitude are required. 

Lower bending stresses, however, are possible for thicker 
viscoelastic plates [ DeBremaecker , 1977], and a regional 
compressive stress has been shown to be unnecessary for 
several forms of plate rheology [e.g.. Chappie and Forsyth, 1979] . 
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It is likely that boundary forces at' zones of active 
continental convergence such as the Alpine-Himalayan belt are 
not of the same form as along oceanic subduction zones. The 
absence of a substantial length of slab and the elevated 
topography at the collision zone suggest that such boundaries 
may offer a net resistance to further convergence. Detailed 
models of the Himalayan zone [ Bird , 1976]give 1 kbar net 
compressive stress due to the load of the mountain belt. 


Drag Forces 

Viscous shear forces will be present at the base of the 
plates if the lithosphere is in motion with respect to the 
asthenosphere . The magnitude of such forces will depend 
upon plate velocity and asthenosphere viscosity and the 
directions will depend upon the velocity of local mantle flow 
relative to the plate. The mantle may be convecting and 

dragging the plates along the tops of convection cells; it 
maybe passively resisting the motion of the plates; or it 
may have a complex flow pattern perhaps only loosely related 
to current plate geometries .and surf ace. motions . These 
possibilities are intimately tied to the questions of the 
radial extent of upper mantle convection and to the depth 
dependence of viscosity in the Earth. 

Convection in the earth was suggested as early as Holmes 
[1929] as the driving mechanism for continental drift. The 
simple early notion that the plates were dragged along the 
tops of large Rayleigh-Benard cells can be shown to be inadequate 
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for plates of large horizontal extent [ Richter , 1973] , though 

this result depends on the assumed depth extent of convection. 

Mantle 'plumes,' are another form of convection proposed as 

a plate driving mechanism [ Morgan , 1972] , though the identification 

of 'hot spots' with plume traces leads to geometries not 

compatible with plumes playing a major driving force role for 

most plates. A major objection to sublithospheric convective 

shear as a major driving mechanism is the amount of time 

7 9 

required for the flow pattern to change, 10 to 10 years 
[Richter, 197 3] , compared to times of 10^ years or less over 

which substantial changes in spreading rates and directions occur 
[e.g., Sclater and Fisher , 1974] 

If the lithosphere is moving with respect to a generally 
passive asthenosphere , there will be a viscous force resisting 
plate motion. The force will be in the direction opposite to 
absolute plate motion and the magnitude of the drag will depend 
on the plate speed and on the law relating shear stress to strain 
rate. .. The simplest assumption is that viscosity is independent 
of strain rate, so that the drag force is linearly proportional 
to the plate velocity. In the actual mantle, a non-linear 
relationship between stress and strain-rate seems likely 
[ Stocker and Ashby , 1973], and -shear heating may also be 
important [ Schubert et al . , 1976] . 

The asthenosphere viscosity may also vary spatially, with 
potentially higher viscosity beneath continental lithosphere 
compared to oceans [ Jordan , 1975] or beneath old oceanic litho- 
sphere compared to young seafloor [ Schubert et a l., 1976]. 
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Most likely the asthenosphere is neither passive nor 
actively dragging the surface plates along their present 
directions of movement, but rather has a complicated flow 
pattern with respect to the surface plate geometries. On a 
large scale, there must be flow to balance the creation and 
destruction of plates. There are several solutions to this 
asthenospheric ’ counterf low' [ Harper , 1978; Chase , 1979; Hager 
and O'Connell , 1979] , depending principally on the assumed 
penetration depth of upper mantle convection. The direction 
of the counterflow is, for at least a few plates, substantially 
different from antiparallel to the 'absolute' plate motion vector 
[e.g., Chase , 1979; Hager and O'Connell , 1979]. An alternative 

approach to the problem of large scale drag is to find that drag for 
each plate that balances the torques due to boundary forces [ Davies , 
1978]. There may also be a smaller secondary scale of 
convection [ Richter and Parsons , 1975; McKenzie and Weiss , 

1975] associated with mantle heat transfer, though this result, 
too, is sensitive to the question of the radial extent of 
upper mantle convection. 

Thus drag forces are among the most difficult to specify 
with any confidence. Shear stress associated with secondary scales 
of convection that may not be important for. net momentum 
transfer to large surface plates may. still be important for 
intraplate stress. At best, the calculation of intraplate 
stress can serve as a test of any specific mantle flow model. 

As knowledge of the flow pattern improves, this test can 


be refined. 
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Transform Fault Forces . 

Shear stress on transform faults separating two plates 
will resist the relative motion of the plates. The seismic 
activity observed on transform faults is a clear indication 
that transform faults resist plate motion. The magnitude 
of that resistance, however, is less clear. If shear stresses 
in the kilobar range exist along major transforms and thrust 
faults at subduction zones, the forces of ridge pushing and 
slab pulling are insufficient to drive the plates and a driving 
mechanism involving shear by mantle convection must be involved 
[ Hanks , 1977] . 

Laboratory measurement o f the shear strengh of rocks at 
pressures and temperatures corresponding to crustal depths 
indicates that several kilobars of shear stress are required 
to produce failure of samples, even when pre-existing faults 
exist [ Brace and Byerlee , 1970; Stesky et al ., 1974]. The 
requirement of high shear stresses is relaxed if the ambient 
effective stress is lowered by high fluid pore pressures. 

Shear stresses in the hundreds of bars range on transform 
faults are indicated, however, by several lines of reasoning. 

Heat flow measurements across the San Andreas fault [ Brune et al . , 
1969], along with' a model of shear heating, indicate a maximum 
shear stress of a few hundred bars acting on the fault. The 
generally orthogonal relationship between transform faults and 
spreading centers has been used to conclude that resistance to 
motion along transforms is less than the resistance to plate 
separation at spreading centers [ Lachenbruch and Thompso n, 1972]; 
to the extent that ridges are the site of upwelling of material of low 
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strength, this result implies that the shear stresses on transforms 
are small. From the inversion of observed plate motions, 

Forsyth and Uyeda [1975] found that shear forces on transform 
faults were not an important part of the driving mechanism, 
because plate velocities do not depend on the length of trans- 
form for each plate. For these reasons and because the total 
length of transform boundaries is less than ridge or trench 
boundaries, we shall adopt the view here that transforms are a 
relatively minor component of the driving mechanism. 


Other Sources of Stress 

Central to the idea that the intraplate stress field can 
be used to study the driving mechanism is the assumption that 
the effect of sources of stress other than plate driving forces 
can be removed or neglected. Such other sources of intraplate 
stress include thermal stresses associated with a cooling 
lithosphere [ Turcotte and Oxburgh , 1973] , stresses arising from the 
motion of an elastic plate on an ellipsoidal earth I Turcotte 
and Oxburgh , 1973] , inhomogeneities in crustal thickness [ Artyushkov , 
1973], lithospheric loading by glaciers, sea level changes, 
and volcanic constructs [ Walcott , 1970; Watts and Cochran , 1974], 
sedimentation and erosion [Vcight and St. Pierre, 1974; Haxby and 
Turcotte, 1976], and residual strains [Swolf s et al . , 1974; Tullis , 


1977] . 




Tensile thermal stresses parallel to the ridge may arise 
as the plate spreads away from the ridge and cools. This stress 
has been estimated to be potentially as large as 40 kbar in 
magnitude [ Turcotte and Oxburgh , 1973] . Thermal stresses, however, 
must be relaxed by deformation on the grain size .level because 
of anisotropic thermal expansion of mineral grains, a process 
that should also relax the larger scale thermal stress without 
coherent lithospheric failure and large intraplate earthquakes 
[ Solomon et al . , 1975] . 

As a plate moves in a north-south direction, it will 

m 

experience a change in radius of curvature due to the ellipticity 
of the earth. This effect will be at a maximum at a latitude 
of ±45° , where estimates of the magnitude of ensuing tensile 
stresses at the edge of the plate, based on thin shell 

theory, are as large as 6 kbar [ Turcotte and Oxburgh , 1973] . 

Most of the earth's plates have velocities with large east- 
west components compared to north-south motions, which would 
minimize the effect of ellipticity. For northward plate 
velocities of a few centimeters a year, it would take 
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about 10 years for the plate to travel from the equator to 45°N. 

In such a time, the plate may deform plastically and relieve 

these membrane stresses. Such stresses may have been important 
in Africa, however; the East African rift system may have been 

generated by the slow northward motion of the African 

plate over the last 100 million years [ Oxburgh and Turcotte , 

1974] . 

Horizontal variations in the density of the crust or mantle 
will result in stresses in the lithosphere. The stresses gene- 
rated by the topographic relief of mid-ocean ridges and continental 


convergence zones arise in this way. _ Features such as mid- 
plate mountains, large plateaus, and continental shelves 
also produce stress, even if the structures are isostatically 
compensated [ Artyushkov , 1973] . Intraplate stress observations 
as indicators of plate driving forces should be well removed 
from such features. 

The addition and removal of mass from the surface of the 
lithosphere will also give rise to stress. Important processes 
capable of producing measurable stress in this fashion are 
glaciation, addition of volcanic constructs, sedimentation and 
erosion. For the case of erosion, the thermal stress associated 
with material cooling as it approaches the earth's surface is 
opposite in sign and may be comparable to the unloading stress 
[ Haxby and Turcotte , 1976] . Evidence on the driving mechanism 
from intraplate stresses should be sought from sites where none 
of these processes are likely to dominate the stress field. 

There is evidence that strains can be locked into rocks 
after the forces producing the strains have ceased to operate 
[ Swolfs et al ., 1974; Tullis , 1977]. Known as residual strains, 
these may be the result of ancient tectonic processes and may 
mask strains associated with present tectonic forces. The magni- 
tudes of residual - strains can be estimated using double overcoring 
techniques. Uncertainties in the origin of the stress field may 
thus be reduced if residual strains are found to be small. Further, 
if the intraplate stress field is uniform over a large area including 
regions of very differeng geologic history, the effects of residual 
strains are probably small. The possibility of some regionally 
consistent remnant stresses indicative perhaps of an old plate 


tectonic event cannot, however, be dismissed. 

EVIDENCE FROM PLATE VELOCITIES 

A primary data base for testing the relative importance 
of the various forces acting on the edges and the bases of 
the plates has been the set of plate motions . Information 
on these forces has been derived from (i) the prediction of 
absolute plate velocities given the relative motions, 

(ii) statistical correlations between plate speeds and force 
distributions, (iii) the prediction and/or the inversion of 
relative plate velocities, and (iv) several of the above at 
geologic times different from the present. 

The motions of the plates with respect to the underlying 
often termed their 'absolute' motions, are derivable 
from plate driving force models as long as the set of forces 
includes some that are dependent on absolute plate' velocity . 

The latitudinal component of predicted absolute velocities 
may be compared with paleomagnetic [ McElhinny , 1973; Jurdy and 
Van der Voo , 1974] or paleosedimentation [ Cla g ue and J arrard, 
1973; van Andel , 1974] data, and the full velocity vector 
may be compared with the motion predicted from the fixed hot- 
spot hypothesis [ Morgan, 1972]. Global calculations of abso- 
lute motions have the advantage that plate boundary forces 
symmetric about the boundary do not enter into the net torque 
balance for the lithosphere [Solomon and Sleep, 1974] . 

V 

Global models of absolute plate models have been determined 
assuming a variety of relationships for li thosphere-astheno- 
sphere interactions,- including a uniform viscous drag lav; 
[ Solomon and Sleep , 1974; Lliboutry . 1974]; enhanced drag 



beneath continental lithosphere; viscous drag balancing slab 
pull concentrated on the subducted plate at convergence zones 
[ Solomon and Sleep , 1974; Solomon et al . , .1975]; minimization 
of the horizontal translation of ridges and slabs [ Kaula , 1975] ; 
and viscous drag characterized by a non-linear relationship 
between stress and strain rate [ Solomon et al . , 1977b] . The 
models with no net torques due to drag also correspond to a 
minimization of the energy dissipated by drag [ Solomon et al . , 

1975]'. For virtually all of these absolute motion models, 
the predicted velocities are indistinguishable on the basis 
of present data and are similar to the velocities calculated 
in a fixed hot-spot frame [ Minster et al . , 1974] . Thus global 
models of present absolute plate motions do not yield much 
information about driving forces except to provide alternatives to 
deep mantle plumes as a physical basis for a hot-spot reference frame. 

Present root-mean-square absolute speeds of the plates show severa 
interesting correlations with boundary or force distributions 
that may have physical significance. The plate speeds generally 
decrease with increasing area of continental lithosphere 
[ Minster et al ., 1974] and increase with the percentage of 
plate boundary connected to a subducted slab [ Forsyth and Uyeda , 

1975] . The plate speed has also been recently shown to be 
generally proportional to the sum of the percentages of plate 
boundary occupied by ridges and connected to downgoing slabs 
[ Aggarwal , 1978] . 

Direct relationships between the relative plate motions 
and driving force models have been sought both as forward and 
inverse problems. A fair fit to observed motions was obtained 
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by Harper [1975] using a ridge pushing force, subduction zone 
pull, and resistive drag; and an improved fit if subduction 
zone forces depend on the age of subducted lithosphere. Formal 
inversions of relative angular velocities to obtain the magni- 
tudes of parameterized driving forces have been presented by 
Forsyth and Uyeda [1975], Chappie and Tullis [1977] and Harper 
[1978] . The solutions of both Forsyth and Uyeda [1975] and 
Chappie and Tullis [1977], using slightly different formulations, 
indicated that the negative buoyancy of the downgoing slab 
and the resistance to relative motion of slab and mantle 
dominate the force balance, with drag beneath continental 
lithosphere and pull toward trenches on overthrusting plates 
small but significant. Harper [1978] , with a substantially 
different formulation . of the forward problem including astheno- 
spheric counterflow, obtained a solution in which subduction 
zone pull, weighted for convergence rate and subducted plate 
age, and continental 'push' due to differential heating are 
the major significant driving forces. 

These same relationships between absolute or relative 
motions and driving force models can be tested for the plate 
system in geologic times different from the present. Global 
absolute motion models for the early Tertiary [ Solomon et al . , 
1977a; Jurdy , 1978] suggest that the negative correlation 
at present between plate speed and continental area and the 
positive correlation between plate speed and percentage of 
subducted boundary may not be generally valid. That dominantly 
continental plates do not always move slowly as at present has 
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also been demonstrated from the rapid pre-Cenozoic velocities 
of several continental masses with respect to the paleomagnetic 
pole [ Gordon et al . , 1978b] . Whether hot spots can be demon- 
strated to define a stable reference frame over periods of 
tens of millions of years depends critically on past plate 
reconstructions [Burke et al ., 1973; Molnar and Prancheteau , 
1975; Solomon et al ,/ 1977a]; absolute motions for the Pacific 
plate derived from a driving force model with ridge and trench 
forces and resistive drag are consistent with a fixed Hawaiian 
hot spot [ Gordon et al . , 1978a] . 

In summary, the efforts at testing models for plate forces 
against relative and absolute plate motions lend confidence to 
the abstraction of the driving mechanism into boundary forces 
and basal drag stress, but leave without definitive answers 
the relative importance of most of the forces in question. 
Testing these same types of force models against the intraplate 
stress field should help narrow the suite. of possible driving 
forces. 


I 



36 . 


TECHNIQUES FOR DETERMINING INTRAPLATE STRESS 

There are a variety of independent methods for inferring 
the local state of stress within a plate, including the 
mechanisms of intraplate earthquakes, direct in situ measure- 
ment, and making use of recent geological features sensitive 
to ambient stress. For most of these techniques, the directions 
of the principal stresses may be recovered more reliably than 
the principal stress magnitudes. The principal stress 
directions form the primary basis for describing 
the intraplate stress field with which calculated stresses 
from driving force models will be compared. 

Earthquake Source Mechanisms 

The source mechanism of an earthquake contains information 
about the orientation of the principal stresses at the source. 

A major advantage of midplate earthquake mechanisms for 
inferring intraplate stress directions is that a medium to 
large earthquake is responding to the average state of stress in 
a volume perhaps several tens of kilometers in dimension. 

If an earthquake is idealized as brittle fracture along a 
fresh, frictionless fault, the maximum (P) and minimum (T) 
compressive principal stress axes will bisect the angles 
between the two nodal planes and will be orthogonal to the 
null axis. Because the seismic wave radiation pattern is 

fundamentally a measure of strain release rather than ambient 
stress and because the assumptions linking the fault orientation 
to principal stress directions are not likely to be generally 
valid, the error in assuming that the principal stress axes 
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bisect the focal planes may be as large as. ±45° [ McKenzie , 1969b] . 

If a large intraplate region is characterized by consistent 
fault plane solutions, the uncertainties may be considerably 
smaller than the maximum of ±45°, however. There are probably 
many pre-existing faults with a distribution of orientations 
in any large area. The inferred principal stress directions 
for earthquakes on pre-existing faults should still reflect 
the ambient stress field for regions with a distribution of 
such faults, since the earthquake will presumably occur 

on that plane which minimizes the resistance to failure. It 
is reasonable, therefore, to assume that the effect of pre-existing 
faults may be neglected for large areas with consistent earth- 
quake mechanisms. 

An isolated earthquake, of course, even if the inferred 

« 

principal stresses are correct, may not reflect the regional 
stress field. The local stress field that caused the earthquake 
may be caused by heterogeneities in topography or loading 
history. 

Most intraplate earthquakes have hypocentral depths of less 
than 20 or 30 km. While 30 km is deeper than most other stress 
measurement techniques can sample, stresses inferred from 
intraplate earthquakes still may not represent the state of 
stress across the entire thickness of the plate. At trenches. 



for instance, the large bending moments associated with the 
downgoing- slab may locally produce marked changes between the 
state of stress at the top and bottom of the plate [ Engdahl 
and Scholz , 1977; Chappie and Forsyth , 1979]. For intraplate 
regions away from topographic loads such as Hawaii, the 
bending moments are probably negligible [ Turcotte , 1974] , and 
it will be assumed that well-characterized fault plane 
solutions represent the stress throughout the plate thickness. 

Fault plane solutions provide information only about the 
orientation, and not the magnitude, of the stresses at 
the source: Other seismic parameters such as apparent stress 

and stress drop may contain some information about the mag- 
nitude of deviatoric stresses, but these quantities are only 
lower bounds. On the basis of available stress drop and 
apparent stress data for intraplate earthquakes, the typical 
shear stress in the lithosphere is at least 100 bars 
[ Kanamori and Anderson , 1975; Richardson and Solomon , 1977] . 

In Situ Methods 

There are two major types of in situ techniques used to_ 
estimate stress. The first, strain relief, uses the relaxation 
of a surface after stress has been removed to calculate the 
original strain [e.g., L eeman , 1971]. With additional infor- 
mation about the elastic properties of the sample, the stress 
can be determined. The second technique, hydrofracture, 
measures stress directly by applying pressure to a closed 
section of a borehole until failure occurs [Haimson and 


Fairhurst, 1967, 1970]. 


The methods for in situ measurement 
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of stress, including their relative advantages and disadvantages, 
have recently been reviewed by McGarr and Gay [1978] . 

Both types of in situ measurements can yield the 
magnitude and direction of the principal stresses. Hydro- 
fracture has advantages over strain relief methods in that elastic 
properties do not have to be measured and in that the technique 
can be applied at great depth. In the Michigan basin, for 
instance, hydrofracture has been used to depths exceeding 5 km 
[Haimson, 1978]. Thus the hydrofracture technique is much less 
sensitive to the local topography that may mask tectonic stress near 
surface. It should be noted, however, that 5 km is still relatively 
shallow compared to the thickness of the lithosphere and such 
measurements need not represent a vertical average of the state 
of stress across the entire plate. One disadvantage of hydro-, 
fracture compared to strain relief methods is that data inter- 
pretation for the former technique is based on the premise that 
one of the principal stress directions is vertical, an assumption 
that may not always hold at depth despite the boundary condition 
at the free surface. 

Interpretations of both strain relief and hydrofracture 
data require considerable care. Stresses due to local topography 
and residual stresses preserved in a medium after the original 
stress has been removed may mask regional stresses, particularly 
for the strain relief data which . typically are taken near the 
Earth's surface or in mineshafts. Another source of uncertainty 
arises when both of the horizontal principal stresses are 
approximately equal. In this case, the azimuth of either stress 
axis is not well constrained. 
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Geological Information 

There are a variety of geological features, ranging 
in dimensions from mineralogical grains to large structures 
on the tens of kilometers scale, which are strongly con- 
trolled by deviatoric stress and which therefore provide 
information on either the orientation or the magnitude of 
stresses acting at the time that the feature formed. 

Information on stress orientation may be obtained from 
the strikes of large-scale faults, folds and rifts [e.g., 
Ahorner, 1975; Johnson and Page , 1976] , from the orientation 
of dikes and of flank volcanic eruption features [e.g., 
Pollard and Muller , 1976; Nakamura et al . , 1977] , and from 
the orientation of joints, cracks, and stylolites [e. g . , 
Greiner , 1978; Kohlbeck and Scheidegger , 1977]. 

Linear island chains have also been interpreted as 
tensional cracks extending through the lithosphere [ Turcotte 
and Qxburgh , 1973] , and thus indicative of principal 
horizontal stresses on a regional scale. Island chains 
with clear age progressions along the length of the chain, 
such as the Hawaiian -Emperor chain, have more commonly 
been regarded as the volcanic products of "hot spots" 
more or less fixed with respect to the deep mantle [ Wilson , 
1963; Morgan, 1972]. Some island chains, such as the 
Line Islands in the Pacific, however, do not show an age 





progression along the length of the chain 
Jarrard and Claaue , 1977], so that a sing 
seems an unlikely explanation. Regardless 
Islands, they are at least SO million 
not likely to be closely related to the c 
pattern. 

Information on the magnitudes of shea 
lithospheric mantle have recently been ot 
microstructure of mantle-derived xenolith 
from the density of dislocations in olivi 
[Goetze and Kohls tedt, 1973] and from the 
xenolith subgrains and recrystalized grai 
1977; Mercier, 1978]. Used in connectior 
geobarometer, such techniques offer the £ 
profiles of differential stress versus d« 
differences so estimated range from a fe*. 
hundreds of bars, though the time period 
stresses are representative and the rela 
stresses to the eruption process that 
to the surface are not well established. 

There are clearly many uncertainties 
of stress measurements. Measurements h< 
relevance for regional stresses when th< 
large areas and where they are corrobor 
independent techniques. In such cases, 
sources of stress is probably small. F 
are several large regions on the earth 
techniques present a consistent picture 



principal stresses in the lithosphere. The data for most 
of these regions, described in detail in the next section, 
are consistent with the inference that these stresses are 
primarily the result of plate tectonic forces. 



REGIONAL DESCRIPTION OF THE INTRAPLATE STRESS FIELD 

Using the above techniques for inferring the directions 
of principal stresses in the lithosphere, we may summarize 
the observed intraplate stress field by region. The data 
used for this summary are listed in Table 1. 

African Plate 

Eastern and southern Africa are two regions where 
the state of stress in the lithosphere has been well studied. 
Geologic evidence, earthquake mechanisms, and in-situ measure 
ments suggest different, but consistent patterns for the 
stresses in these two areas (Figure 1) . 

The tectonics of East Africa are dominated by the east 
African rift system (Figure 1) , which extends from a few 
degrees north to at least 15° or 20 °S [ Fairhead and Girdler , 
1971; Scholz et al . , 1976] . Normal faulting mechanisms for 
earthquakes in the area confirm that the rift valley is an 
active extensional feature, with least compressive principal 
stress trending approximately E-W [ Sykes , 1967; Maasha and 
Molnar , 1972] . One earthquake located near the eastern 
section of the rift southeast of Lake Victoria at 4°S, 35 °E 
shows mainly strike-slip faulting with maximum compression 
trending nearly E-W [ Fairhead and Girdler , 1971] . This 
earthquake, the only known strike-slip event associated with 
the rift system, is located near a normal fault [Fairhead 


and Girdler, 1971] , and would be consistent with motion on a 
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left lateral transform fault between rift centers. (Note that 
Maasha and Molnar [1972] give a different, predominantly 
thrust, mechanism for this event, however.) 

. The western branch of the northern rift system appears 
to be much older, with faulting and grabens that are at 
least Jurassic in age. [ McConnell , 1967]. This branch 
clearly predates the opening of the Red Sea and 
the Gulf of Aden, suggesting that at least the western 
branch is not a newly forming plate boundary [ Girdler and 
Styles , 1974, 1978]. The western branch is still active, 
however, with normal faulting earthquakes indicating minimum 
compressive principal stress trending nearly E-W [ Sykes , 

1967; Banrjhar and Sykes , 1969]. 

In the southern section of the rift system, the 
eastern and western branches merge and there is a general 
decrease in the extent of extension , with the rift valley 
becoming narrower and the walls lower. Normal faulting 
earthquakes, consistent with extension trending E-W, occur 
both within the rift and to the west in Zambia [ Sykes , 

1967; Fairhead and Girdler , 1971; Maasha and Molnar , 1972] , 
though at least' two of the 'largest of these events may be 
associated with the loading by Lake Kariba [ Sykes , 1967] . 

By 20°S, the surface expression of the' rift system dies 
out. A microearthquake study in Botswana between 18 0 and 
21° S indicates that the rift system may be in the early 
stages of propagating further southward [Scholz et al. , 


1976] . 
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While the rift system is clearly the product of regional 
extension, the question arises as to whether or not the 
rift should be classified as an intraplate phenomenon. If 
the rift system is the continental expression of a spread- 
ing center extending from the Afar triangle southward into 
Africa, then the state of stress in the region may not be 
related to forces acting on distant plate boundaries. 

Chas e [1978a] and others have successfully interpreted 
differences in spreading rates between the Red Sea and 
the Gulf of Aden in terms of separation of African and 
Somalian plates along the East African rift system. Inter- 
pretation of the extensional nature of the rift system in 
terms of intraplate stresses is therefore questionable. 

The only large earthquake in eastern Africa located 
well off the rift system has a strike-slip mechanism [Sykes 
and Sbar , 1974] with the P-axis trending NW-SE (Figure 1) . 

The orientation of the P-axis implies horizontal compression 
roughly normal to the trend of the rift system and suggests 
that the state of stress may be very different in areas 
of east Africa away from the rift. 

Southern Africa has a different tectonic setting than 
east Africa. - South of 20°S observations of stress indicate that 
maximum compression trends E— W to NW— SE. These observations 
include earthquake mechanisms and in-situ measurements. 

An earthquake on 29 September 1969 located at 33° S, 

19 °E, near Ceres, South Africa, had a strike-slip mechanism 
with maximum compressive principal stress inferred to be 
nearly E-W [Fairhead and Girdler, 1971; Haasha and Molnar , 


1972]. Aftershocks for this event trended WNW-ESE in 
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agreement with the fault plane orientation. 

In-situ measurements of stress in southern Africa have 
been made primarily using strain^relief methods, generally 
in mines [ Gay , 1975, 1977], As noted above, for such 
measurements it is not always possible to remove the effects 
of ore veins and of nearby cavities and tunnels. Nonethe- 
less, the measured principal stresses show generally consistent 
orientations among measurements at a single site and between 
nearby sites [Gay, 1977] . The directions of maximum horizontal 
compression at six sites listed by Gay [1977] are shown in 
Figure 1. The magnitudes of the deviatoric stresses for 
these sites are up to several hundred bars [ Gay , 1977]. 

In the western part of the African plate, the mechanisms 
of two oceanic intraplate earthquakes [ Sykes and Sbar , 

1974; Richardson and Solomon , 1977] and an in-situ strain- 
relief measurement [Hast, 1969] indicate a NW-SE trend for 
the maximum compressive stress (Figure 1). Sykes [1978] gives 
a thrust mechanism for the earthquake of 23 September 1974 in 
west Africa; the P axis trends roughly NE-SW, but is not well 
constrained . 

Indian Plate 

The Indian plate and relevant continental boundaries are 
shown in Figure 2. Stress directions inferred from 




focal mechanisms for six earthquakes are shown for contin- 
ental India. Five of these solutions are from Chandra [1977] 
and give dominantly thrust faulting; strike-slip mechanisms 
are possible for two of the events, but the inferred P axis 
is not a strong function of the uncertainties in the nodal 
plane solutions. On 10 December 1967 a large earthquake 
occurred near the Koyna Dam at 17.5°N, 73.8 °E. The inferred 
P-axis trends NNW-SSE [ Singh et al ., 1975; Langston , 1976]. 
This event occurred in a relatively aseismic area of the 
Deccan Traps and the P-axis is fairly consistent with 
other events in India. Thus, although the earthquake 
may have resulted from the filling of the reservoir, 
perhaps through an increase in pore -pressure and a reduction 
in effective stress, 'the stress directions inferred from 
the fault plane solution may still reflect the prevailing 
regional stress field. 

In the Indian Ocean basin, much of the seismicity far 
from plate boundaries occurs in the general vicinity of 
the Ninetyeast . ridge, probably an ancient transform 
fault [ Sclater and Fisher , 1974] . Two recent large 
events near the ridge (Figure 2), on 25 May 1964 and 10 
October 1970, have strike-slip focal mechanisms with 
the P-axis trending NW-S E to NNW-SSE, respectively 
[ Sykes , 1970; Fitch , 1972 ] . Thrust fault earthquakes 
occurred to the west of the ridge, on 25 June 1974, and 
to the east of the ridge, on 26 June 1971. The P-axis 
for both events trends NW-SE [Sykes and Sbar, 1974; Stein 
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and Okal, 1978] . The style of faulting in the Indian Ocean 
appears to vary from strike-slip along the Ninetyeast 
ridge to generally thrust on either side of the ridge 
[ Stein and Okal , 1978]'. The orientation of the P-axis 
for all of the events, however, trends approximately NW-SE. 

The seismogenic intraplate stress for all these events 

may thus be similar, with the variation in style arising from 

the availability of a zone of ' lithospheric weakness, the 

Ninetyeast ridge, along which failure may occur. Since 

the Ninetyeast ridge is a major topographic high, the 

stress due to the topography, even though generally compensated 

[ Bowin , 1973] , may be an important contributor to the regional 

stress field . That the inferred axis of greatest compressive stress 

is not orthogonal to the ridge, however, is a strong argument 

against topographic loads dominating other sources of 

regional stress. 

A variety of stress sensitive measurements indicate 
a consistent, but different, stress field for Australia. 

Fault plane solutions exist for five earthquakes in Aus- 
tralia between 1968 and 1973 [ Fitch et al . , 1973; Stewart 
and Denham , 1974;' Hills and Fitch , 1977]. All events have 
either thrust or strike-slip mechanisms with the direction 
of maximum compressive horizontal stress nearly E-W in 
southern Australia and ENE-WSW to N-S in northern Australia 
(Figure 2). The P-axes all have very small plunge angles, 
indicating that the maximum compressive principal stress 
is nearly horizontal. 

Strain-relief in-situ measurements of stress have been 
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made in Australia and Tasmania [ Stephenson and Murray , 1970; 
Mathews and Edwards , 1969; Endersbee and Hof to , 1963]. The 
maximum compressive stresses average about 200 bars, with 
directions that strike E-W in the north and more NW-SE in 
the south (Figure 2), and are only roughly consistent with the 
general trends indicated by the fault plane solutions. 

Most of the Indian plate far from plate boundaries is 
characterized by a state of compressive stress. On the 
Indian continent, south of the Himalayas, the axis of 
maximum compression is nearly horizontal and trends basically 
N-S. The Indian Ocean basin appears to be consistent with 
a more NW-SE trending maximum compression axis and may 
represent a transition from the nearly N-S trend in India 
to the nearly E-W to ENE-WSW trending axis of maximum 
compression inferred from most fault plane solutions and 
in-situ measurements in Australia. 

Eurasian plate 

The state of stress in western Europe has been well 
studied. There are numerous fault plane solutions, in 
situ strain measurements, and a variety of geological indicators 
of past and recent states of stress. A summary of the data 
for western Europe is shown in Figure 3. 
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Fault plane solutions have been determined for a 
large number of primarily strike-slip earthquakes in 
western Europe and have been summarized by Ahorner [1975]. 

Principal stress orientations inferred from these focal 
mechanism data, restricted to earthquakes with published 
fault plane solutions, are plotted in Figure 3. The 
directions for the maximum compressive stress, with few 
exceptions, trend NW-SE. The earthquakes represented include 

events between 1933 and 1971 over a region nearly 500 

* 

kilometers on a side, incorporating parts of the Rhine graben 
and of the Alps. The consistency of the directions argues 
for a large scale source of the stress. Most of the incon- 
sistent data are either pre-1964, and hence predate the World 
Wide Standard Seismograph Network, or are located in the 
geologically complicated Alps. 

In-situ measurements of stress have been made in 
western Europe by many workers and have recently been summ- 
arized by Ranalli and Chandler [ 1975] and Greiner and 
lilies [ 1977] . The magnitudes of the non-lithostatic 
stresses are generally on the order of tens of bars. These 
measurements are typically made using a strain-relief 
technique either near the earth's surface or in tunnels 
or mines. Although such measurements, taken singly, are 
suspect for reasons already mentioned, the consistency of 
the orientations of the maximum compressive stress as shown 
in Figure 3 (open circles) and the agreement with the 
inferred stress directions from mechanism studies are 
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remarkable. 

Other indicators of stress, based on the orientation of 
such geological features as joints, slickensides and 
horizontal stylolites [ Greiner , 1975; Scheidegger , 1977a, b], 
and shown in Figure 3 as open triangles, are also consistent 
with NW-SE maximum compressive principal stress. 

The combination of fault plane solutions, in-situ measure- 
ments, and the orientation of stylolites and major geologic 
features are all consistent with a regional stress pattern 
for western Europe north of the Alps characterized by a 
compressive principal stress which is horizontal and has a 
NW-SE trend. 

There are now focal mechanisms for many earthquakes in 
Asia. Each data point shown in Figure 4 represents the 
average of at least four closely spaced events with 
consistent mechanisms in various seismically active areas 
of Eurasia [ Molnar et al ., 1973; Tapponnier and Molnar , 1977] . 
Molnar et al . [1973] concluded that an "approximately N-S 
to NE-SW trending principal compressive stress appears to 
be transmitted across a large area north and east of the 
Himalayas". The Baikal rift zone, shown in Figure 4, is 
similar to the East African rift system with evidence of 



extensional tectonics C Artemjev and Artyushkov , 1971] . As 
with the East African rift system, it is difficult to 
determine whether such extensional tectonics represent the 
formation of a new plate boundary or are truly intraplate 
in nature. 

The earthquake mechanism data for the Eurasian plate 
east of Europe are presented for completeness, without 
meaning to imply that the solutions can be interpreted 
as strictly intraplate tectonics. Most of the earthquakes 
occur along major faults defining broad zones of extensive 
deformation and folding between relatively stable blocks 
[ Molnar et al . , 1973] . An approach which is based on treating 
these blocks as micro-plates separated by broad boundary 
zones . is misleading because often such boundaries 

fail to completely enclose the stable blocks or die out 
without reaching a micro-plate triple junction [ Tapponnier 
and Molnar , 1977] . The N-S trending compression to 
the north of the Himalayas has been explained in terms of 
the collision of India with Eurasia [ Tapponnier and Molnar , 
1976, 1977] . Both the Baikal extension and the lack of 
comparable deformation in India can also be fit within 
this theory. The major disruption of the Eurasian plate 
by the collision with India suggests that the stresses 
associated with this collision process have been generally 
larger than is typical of intraplate stresses at present. 



Continental Americas 


In South America there have been about a dozen 
intraplate earthquakes for which focal mechanism studies 
have been made. The inferred P-axis for four shallow 
thrust earthquakes and for five shallow strike- 
slip earthquakes located between 200 and 800 km east of 
the South American-Nazca plate boundary [Stauder, 1975] are 
shown in Figure 5. These mechanisms indicate that the 
maximum compressive stress is horizontal and oriented nearly 
E-W at the western edge of the South American plate between 

2 0 and 14 °S. 

Farther eastward in central and eastern Brazil 
there have been several earthquakes that have focal mechanisms 
consistent with NW-SE horizontal maximum compression 
[ Mendiguren and Richter , 1978] . Because these events 
either were small or occurred before the establishment 
of the WWSSN, not all of the focal planes or inferred 
pre-stress axes can be well constrained. Nonetheless, if 
the mechanisms are considered together, a consistent NW-SE 
trend for the P axis is a reasonable interpretation. 

Intraplate stress data are abundant for the North 
American plate (Figure 6) . Perhaps this richness 
contributes to the inferred complexity of the North American 
intraplate stress field. Much of continental North 
America east of the Rocky Mountains and north of the Gulf 
of Mexico is characterized by a horizontal maximum 
compressive principal stress with an E-W to EME-WSW trend 
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[ Sbar and Sykes , 1973, 1977; Haimson , 1977], While there 
are exceptions, the bulk of the evidence', using a 
variety of techniques, strongly suggests a consistent 
regional pattern. The inferred stresses are of probable 
tectonic origin, since local effects of topography and 

glacial history are unlikely to produce such a generally 
uniform trend. It is also unlikely that the compressive stresses 

result from erosion as proposed by Voight- and. St. Pierre [1974] since 
focal mechanism solutions are almost identical for depths from 
1 to about 25 km [ Sbar and Sykes , 1977]. 

The techniques used in North American to infer the state 
of stress at depths from the surface to several kilometers 
include fault plane solutions, hydrofracture and Strain- 
relief in-situ measurements, and such recent geologic 
features as post-glacial buckles or pop-ups. The magnitudes 
of the maximum horizontal compressiye principal stress 
are on the order of a few hundred bars, with values 
varying from about 10 to 600 bars [ Sbar and Sykes ,. 1973] . 

Fault plane solutions exist for a large number of 
earthquakes in North America away from plate boundaries. 

The inferred stress patterns east and west of the Rocky 
Mountains represent different tectonic regimes. North 
America west of. the Rocky Mountains and east of the San 
Andreas fault is characterized by regional extension in 
the Basin and Range Province and normal faulting earthquakes 
along the Intermountain Seismic Belt [ Smith and Sbar , 1974]. 

The origin of these extensional features is uncertain, but 
may well be related to ancient and recent changes in plate 
boundary types along western North America [ Atwater , 1970; 

Burchfiel and Davis, 1975] . The orientation of the T-axis 



for three normal faulting earthquakes along the Intermountain 
Seismic Belt that are typical of the several in the area 
are shown in Figure 7 [ Sbar and Sykes , 1973; Smith and 
Sbar , 1974] . 

Beginning in the Colorado Plateau and extending 
eastward across North America to the Appalachian Mountains, 
fault plane solutions generally indicate thrust or strike-slip 
faulting with maximum compression trending E-W to NE-SW 
(Figure 6) . The data presented in the figure have been 
selected from summaries by Sbar and Sykes [1973] , Raleigh 
[1974] , Hashizume [1977] , and McGarr and Gay [1978] , and 
do not represent all the focal mechanisms that have been 
published. Normal faulting earthquakes on continental shelves 
[ Hashizume , 1973; Sykes and Sbar , 1974] have been excluded because 
the stresses involved in such events may be generated at least in 

part by thermal contraction of the continental margin [ Sleep , 1971] 
or by sediment loading [Walcott, 1972] , and hence may not be related t< 
a regional stress field. Fault plane solutions for four small 
strike-slip earthquakes (Mg = 5. 1-5. 7) in a swarm in the Sverdrup 
Basin in Canada at 77° N, 106° W, not -shown in Figure 6 
because they lie off the map, indicate. NE-SW maximum 
compression [ Hasegawa , 1977] . Additional thrust-faulting 
solutions for North American events north of 60° N were 


reported by Sykes and Sbar [1974] and Hashizume [1974] . 
Finally, whenever several earthquakes in the same locality 
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Although most of the earthquakes indicate E-W to 
NE-SW compression, there are exceptions.' A well studied 
normal fault earthquake occurred on 21 October 1965 in 
southeastern Missouri at 37° N, 91° W [ Mitchell , 1973; 

Street et al . , 1974; Patton , 1976]. The T-axis trends 
NW-SE as shown in Figure 6. The event is located in an 
area of relatively high seismic activity for an intraplate 
region. The mechanisms for small earthquakes in this area show 
changes in fault types and trends of inferred stresses over 
distances of a few hundred kilometers and may well be affected 
by the local structure of the Mississippi Embayment and the Ozark 
Uplift [ Street et al ., 1974; Herrman and Canas , 1978], though 
the better determined mechanisms in the Missouri seismic zone 
indicate mostly E-W to NE-SW horizontal compression [ Herrmann , 
1979]. A similar rapid change in type of faulting is 

seen in the Appalachians, which appear to mark the 

eastern edge of the E-W to NE-SW trend of maximum compressive 

stress for most of eastern North America [ Sbar and Sykes , 

1977] . The eastern Missouri and Appalachian earthquakes 
indicate the complexity of the intraplate stress field 
in North America and the need to consider the effects of 
local geological structures in interpreting focal mechanisms 
in terms of regional tectonic stress. 

Numerous in-situ stress and strain measurements have 
been made in North America. The in-situ data represented 
in Figure 6 have been selected from summaries by Sbar and 
Sykes [1973] , Raleigh [1974] , Ranalli and Chandler [1975] , 


Haims on [1977], and McGarr and Gay [1978]. 
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Post-glacial buckles or pop-ups are geologic evidence 
of the recent state of stress. Such buckles have been 
observed in New York, New England, Ohio, and eastern 
Canada [ Sbar and Sykes , 1973]. Often they are associated 
with quarrying where the lithostatic load has been 
recently removed. Many of the features trend NW-SE, 
indicating an axis of maximum compressive stress 
striking approximately NE-SW [ Sbar and Sykes , 1973]. 

The E-W to NE-SW trend in the maximum compressive 
stress indicated by in-situ measurements and recent geological 
features over a broad region between the Rockies and 
the Appalachians is consistent with the results from 
focal mechanism studies. This corroboration lends 
support to the interpretation of a systematic regional 

stress pattern, presumably of tectonic origin. The fact that 
in-situ stress data are consistent with other measures of stress 

is particularly encouraging because in-situ measurements can 
be made at will in regions of interest , and thus offer the 
potential to greatly extend our knowledge of the intrajjlate stress 
field without awaiting the occurrence of rare and unpredictable 
intraplate earthquakes. 

Oceanic Regions 

For oceanic regions, information about the intraplate 
stress field is limited primarily to earthquake focal 
mechanisms. The number of stress orientation data for 
oceanic regions is thus much smaller than for the contin- 
ents. Further, it is generally impossible to corroborate 



the stress directions inferred from oceanic earthquake 
mechanisms using in-situ measurement techniques. Stress 
directions inferred from a single isolated earthquake 
must always be suspect in their significance for the 
regional stress field, particularly in areas of active 
volcanism or of probable large bending stresses associated 
with the flexural response of the oceanic lithosphere 
to topographic loads. 

There are a number of fault plane solutions for earth- 
quakes in the Pacific plate. Several events are so small 
that well constrained fault plane solutions could not 
be obtained [ Sykes and Sbar , 1974] . The limited first 
motion data indicate thrust faulting, with the P axis 
poorly defined (Figure 7) . Two events on 24 September 
1966 and 28 April 1968 at 12° N, 131° W and 45° N, 175° E, 
respectively, indicate thrust faulting with the maximum 
compressive stress nearly horizontal and trending 
NE-SW [ Sykes and Sbar , 1974] . 

Focal mechanisms have been determined for several 
earthquakes beneath the island of Hawaii, including the 
large events on 23 April 1973 and 29 November 1975 
[ Unger et al ., 1973; Ando , 1976]. The 1973 earthquake had 
a predominantly strike-slip mechanism with the P-axis 
trending NE-SW while the 1975 earthquake mechanism in- 
dicated a low-angle normal fault. Hawaii is not typical 
of a stable interior region of an oceanic plate. The 
stress field may be dominated by large bending stresses 



associated with the load of the Hawaiian chain on the Pacific 
plate [Walcott, 1970; Watts and Cochran , 1974; Rogers and Endo , 
1977] as well as stresses within the volcanic edifices [Fiske 
and Jackson , 1970]. Further, the events may be related to 
stresses associated with magma ascent beneath the active Hawaii 
volcanoes. The mechanisms of Hawaiian earthquakes need bear 
no relationship to the Pacific intraplate stress field. 

Most of the Pacific plate far from plate margins and 
active seamount/island chains is characterized by thrust- 
type earthquakes where focal mechanisms can be well constrained 
Limited data suggest that the maximum compressive stress 
may trend generally NE-SW. 

The Nazca plate has had at least three intraplate 
earthquakes. The focal mechanism for an event on 25 November 
1955 at 17°S, 100°W was obtained by Mendiguren [1971] using P- 
wave first motions and surface wave radiation. The P-axis 
is nearly horizontal and trends E-W (Figure 7) . Two other 
nearby events that occurred in 1944 have mechanisms that 
are compatible with the solution for the 1965 earthquake, 
but the focal planes could not b e further constrained 
[ Mendiguren , 1971] . 

The fault plana solution for an oceanic earthquake that 
occurred in the Antarctic plate at 40°S, 105°W, on 9 May 1971, 
indicates a thrust mechanism [ Forsyth , 1973] with a horizontal 
P-axis trending ESE-WNW 

(Figure 7) . Also shown in Figure 7 are oceanic thrust 
earthquakes in the Atlantic near the mid-Atlantic ridge at 
44 °N, 31 °W, and two near the Caribbean plate at 



approximately 20°N, 60°W [ Sykes and Sbar , 1974]. 

Not included as oceanic intraplate events in our 
compilation are either earthquakes in generally young lithosphe 
[ Sykes and Sbar , 1974; Stein , 1978] or earthquakes 
near trenches presumably related to bending of the 
plate as part of the subduction process [ Hanks , 1971;. 

Chappie and Forsyth , 1979]. Stable oceanic lithosphere 
older than 15 or 20 m.y., is generally characterized by 
a nearly horizontal maximum compressive stress [ Sykes 
and Sbar , 1974], with a direction that appears to have 
a regional consistency (Figure 7). Stein [1979] has 
noted that many large oceanic intraplate earthquakes 
occurred on prominent bathymetric features indicative of 
zones of lithospheric weakness. While such zones of 
weakness may result in a large uncertainty in the specific 
local directions of principal stresses, the sense of 
faulting is still indicative of the general directions of the 
principal stresses, and the consistency of inferred 
stress directions from several events on different 
bathymetric features over a broad area is still strong 
support for a regional stress field of tectonic origin. 

Summary 

A summary of the global intraplate stress orientation 
data is shown in Figure 7. The basic data set shown is 
also documented in Table 1. The data shown represent a 

summary of the more complete regional summaries presented 

in Figures 1-6. Some of the complexity of the intraplate 
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stress field illustrated in those figures cannot be properly 
portrayed in a global view. Only representative data are 
shown for some areas where the data concentration is high 
and the neighboring data are consistent. 

Although the global intraplate stress- field is complicated, 
several large-scale patterns can be seen. Much of stable 
North America is characterized by an E-W to NE-SW trend 
for the maximum compressive stress. South American lithosphere 
beneath the Andes and perhaps farther east in the stable 
interior has horizontal compressive stresses trending E-W to 
NW-SE. Western Europe north of the Alps is characterized by 

m 

a NW-SE trending maximum horizontal compression, while 
Asia has the maximum horizontal compressive stress trending more 
nearly N-S, especially near the Himalayan front. The Indian 
plate has a horizontal maximum compressive stress direction 
that varies from nearly N-S in continental India to roughly 
E-W to NE-SW in Australia. Horizontal stresses are 
variable in Africa, but tend to be compressive with an E-W 
to NW-SE trend in the south. and west and extensional with an 
E-W trend in east Africa. Oceanic lithosphere older than 15 
to 20 m.y. is generally characterized by a compressive state 
of stress. 

The consistency, or lack thereof, of the intraplate 
stress field depends on several factors. Each of the tech- 
niques used for inferring stress has uncertainties that are 
often large or unknown. This is further complicated by the 
fact that measurements are often taken very near the earth's 
surface, or in mines and tunnels where the effects of 
nearby structures may be large. The geologic structure at 
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the measurement site may control the state of stress. 

Thus it is not surprising that the inferred intraplate 
stress field is complicated. Perhaps more surprising is 
the fact that general trends for the orientations of 
principal stresses can often be found that are reasonably 
consistent over regional dimensions. We show below 
that these regionally consistent intraplate stress 
orientations can be quite plausibly explained in terms 
of plate tectonic driving forces. 

NUMERICAL MODELING OF GLOBAL INTRAPLATE STRESS 

To compare the global intraplate stress fields 
predicted by various models of the driving mechanism 
with the observed intraplate stresses summarized above, 
a technique for modeling intraplate stress is necessary. 

Our first approach to this problem was based on the solution 
of a finite-difference approximation to the membrane equations 
for a thin, spherical shell [ Solomon et al . , 1975, 1977b; 
Richardson et al ., 1976], Because of the coarse grid 
used in the finite-difference analysis (10 0 x 10 0 in 
latitude and longitude) , and because of some problems 
with the numerical solution of the equilibrium equations, 
a new finite-element analysis is used in this paper to 
calculate global intraplate stress. The finite-element 
technique has previously been applied to the problem of 
intralithospheric stress in a single plate by Richardson 
[1978] . 



The finite element method used in this paper has 
several advantages over the finite difference technique used 
in our earlier work. Irregular boundaries are easily modeled 
with finite element techniques by varying the size and orien- 
tation of the elements, and nodal forces equivalent to surface 
tractions and line forces acting on each element are 
particularly easy to specify without introducing a minor 
imbalance in the torque acting on the lithosphere. Further, 
the finite element procedure does not require the arbitrary 
adjustments to the grid near the coordinate poles to avoid 
the l/sincf> singularity (where 4> is colatitude) that was an 
artifact of the finite difference coordinate system [ Richardson 
et al ., 1976]. 

Finite Element Technique 

If the lithosphere is approximated by a linear elastic 
shell in the membrane state of stress, the equilibrium condition 
between internal deformation and external loads follows from 

} 

the variational principle of virtual work and leads to the 
standard matrix equation for finite elements [ Bathe and Wilson , 
1976, p. 74] . 

[K] [U] = [F] (1) 

where [U] is the matrix containing the unknown nodal displacements, 
[K] is the stiffness matrix which depends on the geometry and 
elastic properties of the structure, and [F] is the matrix 
containing all loads acting on the structure. 



When the dimensions of the stiffness matrix are very 
large, the Irons wave-front solution technique may be applied 
[ Irons , 1970; Orringer , 1974] . This technique, based on 
Gauss elimination and back substitution, avoids prohibitively 
large in-core computer storage requirements by utilizing 
external storage space. This technique was applied to the 
models of the Nazca plate by Richardson [1978] , although 
the total number of degrees of freedom for those models was 
only about 250. 

The earth's lithosphere is divided into 5246 triangular 
plane stress elements as shown in Figures 8-11, where Figures 
8 and 9 show the mid-latitude regions and Figures 10 and 11 
show the north and south polar caps, respectively. Away 
from plate boundaries, each triangular element has dimensions 
of 5x5x7 degrees. Nodal spacing is approximately three 
degrees along plate boundaries. Twelve plates have been 
included and comparison with the finite difference grid used 
in Richardson et al . [1976] shows a definite improvement in 
the spatial resolution both of predicted stress within the 
plate and of plate boundary definition. Each element has three 
nodes, each of which has two in-plane degrees of freedom. 

These degrees of freedom correspond to the latitudinal and 
longitudinal components of displacement in the global coordinate 
system. There are 2625 nodes associated with the 5246 elements, 
giving a total of 5250 degrees of freedom. 

Specification of Material Properties 

Each element is assumed to have a thickness of 100 km. 

The elastic parameters of the elements away from plate boundaries 
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are assumed to be constant. Values for Young's modulus E 

11 2 

and Poisson's ratio v are taken as 7 x 10 dyne/cm and 0.25, 
respectively. The calculated stresses are independent of the 
specific choice for E as long as the Young's modulus is every- 
where uniform. If the Young's modulus varies spatially, only 
the contrast, and not the actual values, will affect the 
calculated stresses. Variations in plate thickness can be 
modelled by decreasing Young's modulus by a factor proportional 
to the reduced plate thickness. 

Normal and shear stresses are probably transmitted only 

inefficiently across a ridge compared to stable lithosphere 

because of the presence of relatively warm mantle material 

beneath the ridge. To model this effect, the Young's modulus 

10 2 

for ridge elements is assumed to be 3.5 x 10 dyne/cm , or 
a .factor of twenty less than Young's modulus for stable litho- 
sphere. For test cases, the calculated state of stress away 
from plate boundaries was not very sensitive to the choice of 
E for the ridge, as long as the contrast between various 
regions was less than about three orders of magnitude. For 
greater contrasts in E, the various regions were effectively 
isolated and numerical problems arose. 

Convergent plate boundaries may also be regions where stresses 
are inefficiently transmitted from one plate to another. At 
subduction zones, some stress transmitted across the oceanic 
plate may deform the subducted slab rather than act on the 
neighboring plate. Also, inter-arc spreading and island arc 
volcanism indicate the presence of relatively hot material 
near the surface. At continental convergence zones, there may 
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be significant shear heating [ Bird , 1976] . As a first order 
model of the elastic parameters at convergent plate boundaries, 
the Young's modulus is rather arbitrarily chosen equal to 
the value at ridges. 

Transform faults must also act as weak regions because 
strain is concentrated in a narrow zone. This may be particularly 
true for transform faults along the major mid-ocean ridges. 

For these transforms, the lithosphere is usually young and 
therefore warmer than old lithosphere. In all of the models 
considered in this chapter, transform faults along the major 
mid-ocean ridges are assumed to have a Young's modulus 
equal to the value chosen for ridge elements. Transform faults in 
continental lithosphere, however, may be nearly as strong as stable 
plate regions. Deformation along the San Andreas fault appears 
to be strongly concentrated within a few tens of km of the 
fault [ Savage and Burford , 1970] . Transform faults which 
are not dilated or compressed may behave similarly to plate 
interiors with respect to normal stress. Shear stresses clearly 
deform transform faults which therefore are weaker than plate 
interiors for such stresses. In a numerical formulation which 
does not include non-linear failure, this weakness might be modelled 
by assuming that the elastic constants are anisotropic along 
continental transform faults such that the values for Young's 
modulus and shear modulus parallel to the strike of the fault 
are less than the values perpendicular to the strike; future 
modeling should include some form of anisotropic elements. 

For models considered in this paper, continental transform 
fault boundary elements are assumed to have the same uniform 



elastic constants as interior elements. 

The plate boundary elements assumed to have a lower Young's 
modulus than typical plate interiors . are illustrated in Figure 12 

Specification of Driving and Resistive Forces 

The parameterization of various possible plate driving 
and resisting forces are as follows (see Figure 13): forces 
at spreading centers and subduction zones are generally assumed 
to act in the direction of relative plate motion and are 
proportional to the length of plate boundary. The relative 
plate velocities are taken from the RM2 model of Minster .and 
Jordan [1978] , except for the Philippine plate, where the 
Eurasian-Philippine angular velocity of Fitch [1972] has been 
used. Although more up-to-date information on the relative 
motion of the Philippine plate is now available [Chase, 

1978a] the results of the global modeling are unlikely to 
change significantly for a change in the orientation of forces 
acting on the small Philippine plate. 

Forces at plate boundaries might more properly be imposed 
perpendicular to the strike of the boundary. For most ridges 
and for some subduction zones, the difference between the 
direction defined by the perpendicular to strike and the 
direction of relative plate motion is small. For computational 
ease, it is much more convenient to specify the direction of 
forces along a plate boundary with a single rotation pole than 
to require knowledge of the strike of the boundary everywhere. 
Future modeling of subduction zone forces might take into 
account the age of subducted lithosphere, the curvature of 



subduction zones, and the length of subducted slab. For the 
present modeling, it seems sufficient to assume that forces 
at subduction zones act in the direction of plate convergence, 
and that the magnitude of such forces may have a dependence 
on the rate of subduction. 

Plate boundary forces are applied at nodal points 
immediately adjacent to and located symmetrically about the 
boundary. For computational reasons which may be justified 
physically, the forces actually applied at nodes on either 
side of the plate boundary act along great circles locally 
tangent to the small circle about the rotation pole rather 
than along the small circles themselves. With such a specifi- 
cation, symmetric forces acting on the boundary will create 
no net torque on the lithosphere. The difference between 
the directions inferred from great and small circle paths is 
less than a few degrees for all distances greater than 20 degrees 
from the pole of rotation. For distances less than 20 degrees, 
the difference increases. Forces defined by great circle 
directions will remain orthogonal to a plate boundary that lies 
on a meridian about the pole of rotation as the boundary 
approaches the pole. This specification should more closely 
approximate the directions of gravitational forces at a ridge, 
for instance, than would forces defined by small circles about 
the rotation pole. 

The driving force F n at ridges is taken always to be symmetric 

I\ 

about the ridge boundary. Ridge forces adopted are shown in 
Figure 14. For most models, the net force F T exerted on the 
surface plates by the subducted slab is assumed to act symmetrical! 
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about the trench boundary (Figure 15). This condition may 
be relaxed, however, by specifying drag forces to balance 
the torque exerted by forces which are concentrated on the 
subducting plate, a feature of some alternative models. The 
magnitude of F^ may vary spatially to include a dependence 
on the rate of subduction. Continental convergence zone forces 
are shown in Figure 16 . 

Drag forces due to motion of the lithosphere with 
respect to the mantle may be specified for each element. The 
drag force F Q is assumed to be proportional to plate area 
and to the velocity of the plate with respect to the 
underlying mesosphere. The direction of the force is opposite 
to the 'absolute' plate velocity v if the force resists plate 
motions , and given by 

F d = -Dv (2) 

where D is a drag coefficient which may vary spatially between 
continental, young oceanic, and old oceanic lithosphere. Drag 
forces are not specified for plate' boundary elements. At plate 
boundaries the motions between the plate and the mantle material 
below are probably very complicated. At least for ridges, 
the mantle material may have a relatively low viscosity and 
hence the drag forces there may be negligible. 

No forces are explicitly applied along transform faults. 

If drag forces dominate the driving mechanism, shear stresses 
along transform faults will be high in the models as a consequence 
of the rapid variation in the direction of drag forces across 
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the transform, and so will be consistent with the inference 
that resistance on transforms is significant for plates driven 
by basal shear [Hanks, 1977]. 
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MODELS OF INTRAPLATE STRESS 

In this section various models of the driving 

mechanism are tested against the observed intraplate stress 

field using finite element analysis to calculate predicted 

stresses. The first class of models includes only 

symmetric forces acting on plate boundaries. In later 

models, drag forces on the base of the plates are 

included. A representative list of all models studied 

is given in Table 2. The various force parameters include 

F R/ the symmetric force exerted at ridges; F^ , the symmetric 

force exerted at trenches; F g , the additional force 

exerted on the subducted plate at trenches; F v , a 

velocity-dependent symmetric force at trenches (see below) ; 

F^, the symmetric force at continental convergence zones; 

and D_, D_, D, , the drag coefficients (equation 2) beneath 
coy 

continental, old oceanic, and young oceanic lithosphere, 
respectively. The forces F R , F T , F g , F v are positive if 
they drive plate motion, and F_, D , D , D are positive 

v w y 

if they resist plate motion. 

Forces only at Plate Boundaries 

The simplest models each include only one of the force 
parameters. While all of these models are too simple to 
match the inferred intraplate stress field, they each 
isolate the effect of one of the possible forces. Because 
we are concerned only with intraplate stress, a trade-off 
for symmetric forces exists between the magnitude of the 
applied forces and the stiffness of boundary elements. 
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For example, additional resistance at a’ trench boundary 
could be included by increasing the stiffness of the boun- 
dary or by decreasing the driving forces. If the strain on 
boundary elements was required to be proportional to the 
spreading (or subduction) rate, this trade-off would no 
longer apply. 

The first such model, El, includes only a horizontal 
pushing force F R , equivalent to 100 bars compressive stress 
across a 100 km thick plate, at all spreading centers 
(Figure 14) . The predicted deviatoric intraplate stress 
field for model El is shown in Figure 17-. Though the 
calculated stresses here and throughout will be refered to as 
deviatoric stresses, they are strictly the differences 
between the horizontal and vertical stresses. In classical 
continuum mechanics, deviatoric stresses are measured 
with respect to the hydrostatic, rather than the vertical, 
stress. We shall nonetheless use the term deviatoric stress 
in this non-standard fashion to emphasize that a calculated 
tensile stress does hot necessarily imply a net tensional 
state of stress in the earth. 

Most of the earth's surface is predicted to be in a 
state of deviatoric compression for model El. The typical 
magnitude of the calculated stresses is about: 50 bars. 

Away from plate boundaries the stress field various smoothly. 
Near ridge boundaries the maximum compressive stress is 
generally aligned with the local direction of plate motion. 
Both of these features are as one would expect, and lend 
confidence to the numerical techniques that have been used. 


73 . 


Comparison with the observed intraplate stress field shown in 
Figure 7 shows moderate agreement in several regions. In 
eastern North America the predicted maximum compressive 
principal stresses trend N-S to NE-SW. In the western portion 
of the Nazca plate a state of deviatoric compressive stress 
is predicted, although the direction is poorly constrained 
because the two principal stresses are of nearly the same 
magnitude. The southern portion of the Indian plate shows 
a NW-SE to nearly E-W trend for the maximum compressive stress. 

The Pacific plate is generally in a state of compression. 

A second simple model, E2, has only symmetric pulling 
forces F t at trenches (Figure 15) equivalent to a deviatoric 
tensile stress of 100 bars across a 100 km thick plate. The 
predicted stress field for model E2 is shown' in Figure 18. 

Most of the earth's surface is in a state of deviatoric tension, 
with typical values of about 50 bars in the middle of the plates. 
The predicted stresses for this model are inconsistent with 
the observed stress field. Although this is a very simple 
model, it indicates that slab pulling forces cannot completely 
dominate the intraplate stress field. 

The net pulling force due to the .slab is the difference 

between the gravitational force and the resisting forces acting 

on the slab as it moves with respect to the mantle. This net 

« 

force may depend on the subduction rate. For subduction rates 
greater than about 9 cm/yr, the gravitational force on the slab 
may be nearly constant [ Richter and McKenzie , 1979] . For lower 
subduction rates there is likely to be a lesser gravitational force 
due to the reduced thermal contrast between slab and mantle. 
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The viscous drag resisting motion of the slab may increase 
with velocity until, at some terminal velocity, the resisting and 
pulling forces approximately balance [ Forsyth and Uyeda , 1975] . 

To model these ideas, a slab force F^ is assumed which increases 
linearly from zero with subduction rate v from v = 0 to 4 cm/yr, 
decreases linearly to zero from v = 4 to 12 cm/yr, and is zero 
at higher v. At v = 4 cm/yr, F^ is equivalent to a deviatoric 
tensile stress of 100 bars across a 100 km thick plate. The . 
predicted stress field for a model (E3) with only F v present 
is quite similar to that for E2, with generally extensional 

m 

horizontal stresses, but with relatively higher stresses near some 
trenches (Tonga-Kermadec , Mariana, Kuril-Japan) . 

As noted above, convergence zone forces are probably 
different for zones of continental convergence such as the 
Himalayan and Alpine belts than for regions of oceanic litho- 
sphere subduction. The force at continental convergence zones 
from the uplifted mountain belts should exert a compressive 
stress on the plates and thus yield a net resistance to 
convergence. As a first approximation, the force per unit 
length along those portions of the African-Eurasian and Indian- 
Eurasian plate boundaries undergoing continental collision are 
assumed to be equal and uniform, in spite of considerable struc- 
tural complexities in both regions. The directions of these 
forces (Figure 16) are nearly perpendicular to the strike of the 
topography associated with the convergence zones. The calculated 
stresses for a model, E4 , with only continental convergence zone 


forces are 



shown in Figure 19, where F c has been assumed to be equivalent 
to a deviatoric compressive stress of 100 bars. The largest 
stresses, about 50 bars, occur in the Eurasian and Indian 
plates near the Himalayas. As expected, the magnitude of the 
stresses decreases away from the plate boundaries. In North 
America, for example, the stresses are only a few bars. 

Comparison between Figures 19 and 7 shows several 
interesting features. First, the directions of the calculated 
maximum compressive stresses in northern India and Asia 
are similar to those inferred from the focal mechanisms for 
these areas. Also, the predicted maximum compressive stress 
varies from N-S to NW-SE between northern and central portions 
of the Indian plate. In the northern portions of the Asian 
plate the component of deviatoric tension is consistent with 
Baikal rift tectonics. The calculated stresses in Asia for 
model E4 are qualitatively similar to the results obtained by 
Tapponnier and Molnar [1976, 1977] from a slip-line analysis 
which treated India as a rigid indentor impinging upon Asia. 
The calculated maximum compressive deviatoric stress in Europe 
for model E4 trends NW-SE, also in agreement with the data 
for this area. 

For models considered later in this chapter, forces 
associated with continental convergence zones appear to be 
important for intraplate stresses in regions adjacent, to such 
zones, but have little effect at great distances. 

For this simple class of models with only plate boundary 
forces, it is not possible to consider forces which are con- 
centrated on one side of subduction zones. Below we consider 



models including a force F<. at trenches which acts only on 
the subducted plate. Such a force exerts a net torque on 
the lithosphere, which must be balanced by drag forces to 
satisfy mechanical equilibrium. 

A linear combination of models is possible because the 
deformation of an elastic lithosphere is linear in the applied 
forces. To consider the effects of both ridge and convergence 
zone forces, the results from models E1-E4 are linearly super- 
posed in models E5 through E16 (Table 2). 

The results of comparing predicted and observed stress 
for these models with only plate boundary forces may be briefly 
summarized. Models that produce a' reasonable- fit to the 
observed data (e.g., E15, E16) all include substantial ridge 
pushing forces. For models with boundary forces only, an 
upper limit of about 3:1 for F./F^ may be set. For higher 
ratios of F /F , the predicted horizontal stresses are dominantly 
extensional and fit the observations poorly (Table 2) . If F v 
(rate-dependent slab forces) replaces F T , the bound F V /F R is 
slightly higher, about 5:1. Forces at continental convergence 
zones bordering Eurasia are important for the intraplate stress 
field in nearby regions, but have insufficient effect in other 
parts of the lithosphere to increase the upper bounds established 
for the ratios F T /F R and F V /F R . 

The boundary forces considered above, with the exception 
of continental convergence forces, are assumed to drive, rather 
than resist, plate motions. Although the net torque exerted 
on the lithosphere is zero for all of the boundary force models 
considered thus far, forces acting to resist plate motions have 


not been explicitly included. The next set of models will 
consider the role of drag forces acting on the base of the 
lithosphere . 

Viscous Drag Forces 

Drag forces on the base of the lithosphere depend on the 
velocity of the plate with respect to the mantle and on the 
mantle viscosity. As discussed above, the motion of the plate 
with respect to the mantle is poorly known. For the first set 
of drag models considered, a linear drag law of the form given 
in equation (2) will be assumed. The absolute velocities are 
determined from the condition that the drag forces exert no 
net torque on the lithosphere [ Solomon and Sleep , 1974] . The 
drag coefficient D is positive if F D resists plate motions. 

Model E17 is characterized by drag forces of the form 
given in equation (2) , where the drag coefficient is taken 
to be everywhere uniform. The calculated absolute plate 
velocity field for the Pacific plate is given by an angular 
velocity of 7.0 x 10 ^ deg/yr about a pole at 62.4°S, 111.1°E, 
in good agreement with other similar models [Solomon et al., 
1975] . Absolute motions for the other plates may be calcu- 
lated from the plate velocities relative to the Pacific plate 
[ Minster and Jordan , 1978] . We should note that model E17 
can never be a physically realizable model by itself, since 
viscous shear at the base of the plates is acting to resist 
motion induced by other forces not explicitly included. We 
present results for models with only drag forces to isolate 
the effects of various forces which are incorporated in later, 


more complete models. 


The calculated stresses for model E17 are shown in Figure 
20. The largest deviatoric stresses are tensile in the 
western Pacific and Asia and compressional in. the eastern 
Pacific and the Nazca plate. The stresses are small in 
eastern North America and Europe. The calculated stresses are 
in poor agreement with the observed stress field in the western 
Pacific, Asia, and the Indian plate. 

The drag coefficient may be enhanced beneath continents 
relative to oceans f Jordan . 1975] or beneath old oceanic 
lithosphere relative to younger regions [ Schubert et al ., 1976]. 
Model E18, with drag only beneath continental lithosphere, and 
model E19, with drag forces only beneath oceanic lithosphere 
older than 80 m.y. , were constructed to explore the effects 

of various parameterizations of the drag coefficient. From 
the comparison of predicted and observed stresses for models 
E17-E19 and their combinations, we may make the following 
conclusions. If only drag forces act on the plates (an 
unreasonable assumption) , the drag coefficient beneath oceanic 
lithosphere must be non-zero if compressive stresses are 
required for oceanic regions. The calculated stresses are 
similar in direction for uniform drag, and old oceanic drag 
force models. The drag coefficient beneath continents may be 
higher than beneath oceans by a factor of up to 10 without 
significantly affecting the continental stresses. 

A model in which drag forces drive rather than resist 
plate motions can be constructed from model E17 by assuming 
that the drag forces act in the direction of absolute plate 
motions. For such a case, the. sense of stress at each point 



in Figure 20 is reversed and the magnitude of the stress is 
unchanged. Such a model . (E20) predicts deviatoric compression 
in Eurasia, the western Pacific and Indian plates and 
deviatoric tension in the eastern Pacific, Nazca and South 
American plates. 

Models with only drag forces cannot be expected to match 
the observed intraplate stresses very well, since the contri- 
bution of plate boundary forces due to ridges, subduction- 
zones, and continental convergence zones have been neglected. 

In the next class of models, both plate boundary and drag 
forces are included with the aims both of establishing limits 
on the various contributions to the driving mechanism and of 
finding a best model. 

Combined Boundary and Drag Forces, 

Models E21-E23 are designed to test whether the limit of 
about 3:1 for F /F n established for the boundary force models 
may be increased by adding resistive drag forces. Stresses 
for model E21, in which a small amount of drag has been added 
to model E7 are very similar to those for model E7, with small 
differences confined primarily to the Pacific plate. The effect 
of increasing the drag coefficient in model E22 is to increase 
the magnitude of the calculated deviatoric tensile stresses 
in Asia and the western Pacific and of the compressive deviatoric 
stresses for the eastern Pacific and Nazca plates. The stresses 
in' North America and Europe are not changed appreciably. 
Calculated stresses for model E23, with a lai-ge drag coefficient, 
include dominantly deviatoric tension in the western Pacific, 
Asia, and southern Indian plates. The predicted stresses are 
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compressive in the Nazca plate, and agree with the earthquake 
data for this plate. The maximum compressive stress is small 
and oriented NW-SE in eastern North America. 

Models E21-E23 indicate that adding uniform drag forces 
does not increase the upper bound of about 3:1 for F T /F R . The 
only areas where adding drag forces improves the fit to the 
data are the western Pacific and Nazca plates. 

Concentrating the drag coefficient by a factor of six 
beneath continents produces stresses very similar to those 
for model E21, where the drag coefficient was assumed to be 
uniform, with the main difference being that the compressive 
component of the stresses is slightly increased for continents 
and slightly decreased for oceanic lithosphere for model E24 
compared to model E21. The fit to the observed data is only 
marginally changed. Even if the drag coefficient beneath 
continents is enhanced by a factor of 11 compared to oceanic 
lithosphere (model E25) , there are only minor changes in the 
stress field (e.g., in South America and the Pacific) compared 
to E24. 

The fit to the observed stresses - in Europe, Asia, and. India 
can be improved over models E21-E25 by including forces at contin- 
ental convergence zones-. (model E26). The fit to the 

observed data is essentially unchanged away from the continental 
convergence zones. 

Models E21-E26 (see Table 2) indicate that resistive drag 
forces acting on the base of the plate cannot significantly 
improve the fit to the observed intraplate stress data for 
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the limiting case of slab to ridge forces. The drag coefficient 
beneath the continents may be larger by a factor of ten 
with little effect as long as drag forces do not dominate 
the intraplate stress field. 

In the previous models, drag forces were assumed to resist 
plate motions. A simple model (E27) in which drag forces 
drive plate motions with resistance to plate motion provided 
by compressive forces at continental convergence zones and 
at trenches predicts stresses as shown in Figure 21 (Table 2) . 

Drag is assumed to be in the direction of local absolute plate 
velocity (D<0) in equation (2). The fit to the data is poor ' 

for the eastern Pacific, Nazca, and South American plates. 

The calculated stresses in eastern North America are acceptable, 
but the orientations of the maximum compressive stresses in " 

Europe and Asia are inconsistent with the data. The calculated 
stresses for model E27, in general, are at best in only 
moderate agreement with the data. 

The effect of increasing the drag coefficient for models in 
which drag drives plate motions (E28) is to predict larger deviatoric 
tension in the eastern Pacific, Nazca, and South American plates. 
Increasing the drag coefficient will further degrade the fit 
between calculated and observed stresses in these regions. 

It is unfair to conclude that viscous drag does not 
drive the plates from simple models such as E27-E28. However, 
for the assumptions made about the forces acting on the plate 
boundaries and about the relationship between drag forces and 
the velocities of plates with respect to the mantle, models 
with resistive drag forces are in better agreement with the 
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observed intraplate stresses than are models with plate-driving 
drag forces. 

It is probably an oversimplification to assume that drag 
forces either uniformly resist or drive plate motions. . .An 
alternative approach, taken by Richardson [1978] and Davies 
[1978] , is to suppose that drag balances the net torque on 
each plate due to boundary forces. The resulting drag will 
vary from plate to plate and will not bear a simple relation- 
ship to relative plate motions in contrast to the drag derived 
from absolute plate motion models consistent with known relative 
velocities. It is straightforward with such an approach to 
consider subduction zone forces that act only on the subducted 
plate. 

A model of driving forces (E29) with symmetric forces 
at ridges and continental convergence zones and with subduction 
zone forces acting only on the subducted plate is considered 
as a test of the Davies [1978] concept. The basal shear or 
viscous drag on each plate arises from the rotation of the plate 
with respect to the underlying mantle (which may be moving) . 

The basal shear necessary to balance the torque on each plate is 

given in Table 3. Note that the poles for drag to balance 
boundary forces do not always correspond closely to hot-spot 
poles [e.g., Minster and Jordan , 1978].' 

Davies [1978] assumed that kilobars of shear stress act 
on transform faults [c.f.. Hanks , 1977] and used the results of 
Chappie and Tullis [1977] to estimate slab forces. Because of 
the large resistive forces assumed at trenches, Davies concluded 



that plates without significant subducted slabs such as South 
America and Eurasia were driven from below by drag forces. 

For model E29, no forces have been specified on transform faults 
or on the non-subducting plate at subduction zones. Thus, for 
example, no drag forces are specified on the Caribbean plate 
Though a direct comparison of the Davies [1978] model 

with model E29 is therefore not strictly possible, model E29 repre- 
sents a first attempt to model subduction-zone forces that act 
only on the subducted plate. 

The stresses for model E29 are shown in Figure 22. For 
several areas the predicted stresses agree very well with the 
data. In the North American and Nazca plates, the orientation 
of the maximum compressive stress is well matched by the model. 

The fit is almost as good in Europe and in Asia north of the 
Himalayas. In the Indian plate, compressive stresses trend 
NW-SE in continental India, in agreement with the data, but 
the fit is poorer in Australia. In South America, the 
maximum compressive stress trends E-W, in only moderate 
agreement with the data. In the Pacific and the eastern 
part of the African plate the agreement with the data is poor. 

On the whole, the model provides a better fit to continental 
than oceanic data, and suggests that any force pulling the 
overthrust plate toward the trench is probably lower in 
magnitude than the net' pull on the subducted plate. 

The last three models represent an attempt to find a best 
model, rather than limiting models, of the driving mechanism. 

Models E30-E32 have symmetric forces at ridges and continental 
convergence zones equivalent to a deviatoric compressive stress 
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of 100 bars (200 bars for continental convergence zones in 
E32) across a 100 km thick lithosphere, symmetric forces at 
subduction zones equivalent to a deviatoric tensile stress of 
100 bars, and drag forces on the base of the lithosphere with 
a drag coefficient that is concentrated by a factor of four 
(E30) to six (E31-E32) beneath continents (Table 2) . The 
predicted directions of principal stresses for model E31 
are shown in Figure 23; stresses for models E30 and E32 are 
generally similar. The model stresses are in good agreement 
with the data for eastern North America, Europe, Asia near 
the Himalayas , and the Indian plate. The fit to the data is 
good in South America, especially away from the trench, and 
in western Africa and is acceptable in most of the Pacific 
plate. ’ The orientation of the calculated maximum compressive 
stress in the Nazca plate for model E31 is only in moderate 
agreement with the orientation inferred from the single fault 
plane solution available [ Mendiguren , 1971] . The fit to the 
data in the northern Pacific, eastern Asia, and east Africa is 
rather poor. The fit to the data in the northern Pacific and 
eastern Asia could probably be improved if subduction zone or 
drag forces were decreased along the western Pacific plate 
margin or if slab forces were concentrated on the subducted 
plate. No attempt, however, has been made to vary plate boundary 
forces locally to match inferred stresses. If such an approach 
were adopted, most observed stresses could probably be matched 
but the solution for the driving mechanism would be unjustifiably 
arbitrary and non-unique. The approach used in this paper has 
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been to treat each potential driving and resisting force acting 
on the plates by class, and make estimates of the importance 
of each type of force. Thus, while models E30-E32 are good 
models in terms of the fit between calculated and observed 
stresses, the actual driving mechanisrrs for plate tectonics 
are not tightly constrained to the particular forces that 
have been assumed for these models. The role of modeling 
intraplate stresses is best suited to testing the relative 
importance of various forces acting on the plates, rather than 
determining the exact balance of forces that make up the 
driving mechanism. 



86 


SUMMARY AND CD NCL US 10 NS 


With a 'variety of techniques it appears to be possible 
to measure reliably the orientations of principal stresses 
within the plates. Over large regions of stable plate interior 
the directions of principal, stresses as inferred from intra- 
plate earthquake mechanisms, in situ measurements, and recent 
stress-sensitive geological features show consistent trends 
indicative of a regional stress field. In areas far from 
recent deglaciation, large-scale topographic loading and 
erosion, or thermal activity, such, regional stresses may be 
reasonably ascribed to the tectonic driving and resistive 
forces associated with the motions of the plates. This 
paper tests models of those forces by calculating the intra- 
plate stress predicted from each model and comparing the 
predicted and observed orientations of principal stresses. 

The calculated. intraplate stress field is very sensitive 
to forces applied at plate boundaries and on the base of 
the plates.’ The predicted stresses are in best agreement with 
the observed intraplate stress orientations when pushing forces 
at ridges are included and when the net pulling force 
of subducted lithosphere is less than a few times larger than 
other forces acting on the plates. Resistive forces associated 
with trench thrust faults and motion of the 
slab with respect to the mantle must nearly balance the 
negative buoyancy of the relatively cool, dense slab [ Smith 
and Toksoz , 1972; Forsyth and Uyeda , 1975] . The upper limit 
on net slab forces may be extended by less than a factor of 
two if slab forces are assumed to depend on the subduction 
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rate such that the force acting on the fastest moving plates 
(e.g., the Pacific and Nazca plates) is reduced. 

Forces resisting further convergence at continental 
collision zones along the Eurasian plate are important for 
intraplate stresses, and improve the fit to the data, in 
Europe, Asia, and the Indian plate. The upper bound on the 
ratio of net slab force to other forces is not affected by 
the inclusion of continental convergence zone forces. 

Resistive viscous drag forces acting on the base of 
the plate improve the fit to the intraplate ' stress field 
for the Nazca and South American plates as long as some drag 
acts beneath oceanic lithosphere. The intraplate stress 
field is not very sensitive to an increased drag coefficient 
beneath old oceanic lithosphere compared to young oceanic 
or continental lithosphere. Increasing the drag coefficient 
beneath continents by a factor of five or ten changes the 
calculated stresses only slightly and has little effect on the 
overall fit of calculated stresses to observed stresses as 
long as some resistive drag acts beneath oceanic lithosphere. 

Models in which drag forces drive rather than resist 
plate motions are in poor agreement with the data. This 
poor agreement may depend on the oversimplified model of 
the adopted interaction between the plate and the 
asthenosphere. The actual flow pattern in the mantle, 
including possible multiple scales of convection [Richter and 
Parsons , 1975] , and the resulting drag forces acting on 

the base of the plates, may be considerably more complicated 
than has been assumed in these models. 


An alternative approach to modeling drag is to require 
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drag forces to balance the torques due to boundary forces 
for each plate [ Davies , 1978; Richardson , 1978]. The drag 
forces thus depend on the assumed boundary forces. The calcu- 
lated stresses for one choice of boundary forces (model E29) , 
with symmetric forces at ridges and continental convergence 
zones and forces exerted only on the subducted plate at 
subduction zones, are in good agreement with the data for the 
Nazca and several continental plates but in poor agreement 
with the data for oceanic plates near subduction zones. 

The calculated stresses for several regions suggest 
limitations of the modeling and indicate improvements that 
might be included in future work. Most models failed to 
predict extensional deviator.ic stresses in east Africa. 

Africa is surrounded by ridge and continental convergence 
plate boundaries which 1 have been assumed in this p^per to 
compress the plates. The absolute velocity of the African 
plate is small compared to most plates, and drag forces have 
been assumed to be proportional to velocity. Increasing the 
drag coefficient sufficiently to affect African intraplate 
stresses degrades the fit to the data in most other regions. 

If the assumed form of boundary and drag forces is correct, then 


some 



89 . 


other forces must dominate the tectonics of east Africa. 

Forces associated with northward motion of Africa on an ellip- 
soidal earth during the past hundred million years may pro- 
duce the extensional tectonics of east Africa [ Qxburgh and 
Turcotte, 1974]. Alternatively, the east African rift may be 
a plate boundary between the African and Somalian plates,, 
rather than an intraplate feature [Chase, 1978a] . Future 
modeling of the intraplate stress field might treat the rift 
system as a spreading plate boundary. 

Some models predict deviatoric tension in the Nazca 
plate while predicting stresses elsewhere in reasonable 

agreement with existing data. It may be that the inactive Galapagos 
Rise, ignored in this modeling, may contribute to the state of compres- 
sion inferred from the thrust earthquakes in the Nazca plate [c.f., 
Richardson , 1978]. Alternatively, the slab forces acting along the Peru 
Chile trench may be less than have been assumed. Support" for 
such a view may be found in the pronounced gap in intraslab 
seismicity between 200 and 500 km depth beneath South 
America [ Isacks and Molnar , 1971] . The portion of the slab 
able to pull on the Nazca plate may thus be confined to depths 
above 200 km, and the net pull may be less than has been 
assumed in the simple models of slab forces used here. 

For most models, the calculated stresses for South 
America and eastern Asia have a large component of deviatoric 
tension that is not supported by the data. This may result 
from the assumption for these models that slab forces act 
symmetrically about the trench axis. A model with slab forces 
only on the subducted plate produced stresses for South 


/ 
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America and eastern Asia in much better agreement with the data 
[c.f., Mendiauren and Richter . 1978], The assumption of symmetric sub 
duction zone forces is undoubtedly an oversimplification and 
methods for balancing the torque due to asymmetric slab 
forces other than simple drag forces acting on each plate 
should be explored. 

The Andes may represent another source for compressive 
stresses in South America. In much the same manner that 
the Himalayas and Alps may stress the plates, the horizontal 
density variation associated with the Andes may produce 
compressive stresses in South America. Future modeling might 
include the effects of major mountain ranges and perhaps the 
effects of horizontal density variations associated with 
continental margins [ Artyushkov , 1973] . 

The number of parameters associated with driving mechanism 
models is already fairly large (e.g.. Table 2). It is likely 
that if even more parameters were considered, such as sub- 
duction zone forces that varied from region to region, the 
fit to the data could be improved. The models are already 
non-unique, and the philosophy behind the modeling in this 
paper has been to describe the dominant features of the 
driving mechanism that are functions of boundary type with 
as few free parameters as possible. 

In future modeling of this kind it would be worthwhile 
to compute plate velocities in addition to intraplate stresses. 
Conceptually, this calculation is straightforward because 
the viscous equation for stress' 
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a . . = Xv . . 6 . . + n (t — — + x — — ) 
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where v^ is -the ith component of velocity and X and n are the 
bulk and shear viscosities, is homologous to the 
elastic equation for stress. Thus, an ordinary finite 
element routine can be used if displacements are interpreted 
as velocities and elastic constants are interpreted as 
the homologous viscous constants. A finite bulk viscosity 
must be assumed, because in two dimensions there is a net 
dilatation or compression of the plates at plate boundaries. With 
a set of viscous constants for different kinds of plate 
boundaries, driving forces would be fit to observed plate 
velocites as well as intraplate stresses. The drag forces, 
if any, then would be automatically self-consistent with 
the plate velocities in a successful model. 

There may be some numerical difficulties with such an 
approach. It would be desirable to impose a large viscosity 
contrast between plate boundaries and plates so that plate 
velocities would be well defined within plates. This contrast, 
however, will cause the velocity to be nearly constant within the 
plates. Round-off errors in computing the horizontal derivatives 
of velocity and hence stresses may thus cause numerical 
problems. A small amount of strain within 

! 

plates will probably have to be tolerated in such models. 

An alternative numerical approach would be to include forces 
due to drag at each element which were proportional to the calcu- 


lated displacement. This approach would require 
cation of standard finite element routines. 
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The main contribution of this work about the question of 
the driving mechanism for plate tectonics is the demonstration 
that observed intraplate stress directions can be explained 
in terms of the forces acting on the edge and along the base 
of the lithospheric plates. Comparison of observed intra- 
plate stresses with calculated stresses for various models of 
the driving mechanism is an effective test of the relative 
importance of plate driving and resisting forces. The modeling 
can be extended in the future to include distributed forces 
due to topographic features such as continents, thermal 
forces associated with regions of anomalously thin or thick 
lithosphere, and explicit . shear resistance along transform 
faults. The role of drag forces can be further tested for 
specific and realizable flow patterns in the mantle. Additional 
data on the intraplate stress field, particularly for oceanic 
regions, should help to define regional patterns that 
are now poorly known. ' Of special importance for further 

constraining models of the plate tectonic driving mechanism 
is the reliable determination of the magnitudes of principal 
lithospheric stresses on a regional scale. 
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11 


3 


ca 39. 

ca- 89.5 

60 


n 

12 


3 


ca 36. 

ca- 98. 

65 


kk 

13 


3 


ca 39.5 

ca- 78.5 

25 


n 

14 


3 


ca 33. 

ca- 80.3 

50 


rr 

15 


3 


ca 30.8 

ca- 98 . 3 

300 


bb 

16 

25Nov65 

1 

5.75 -17.1 

-100.2 

276/0 * 

10/74 

V 

17 

9May71 

1 

6.2 

-39.8 

-104.8 

106/15 

286/75 

j ' 

18 

10May63 

1 


- 2.2 

- 77.6 

250/6 

341/16 

qq 

19 

9Feb67 

1 

6.3 

- 2.9 

- 74.9 

243/15 

340/24 

qq 

20 

- 3Nov63 

1 


- 3.5 

- 77.8 

79/8 

239/81 

qq 

21 

19Jun68 

1 

6.1 

- 5.5 

- 77.2 

272/9 

115/80 

qq 

22 

20Mar72 

1 

6.1 

- 6.8 

- 76.8 

262/20 

76/70 

gg 

23 

14Feb70 

1 

5.8 

- 9.9 

- 75.6 

264/17 

124/68 

gg 

24 

24Jul69 

1 

5.9 

-11.9 

- 75.1 

99/15 

206/48 

gg 

25 

10ct69 

. 1 

5-8 

-11.9 

- 75.1 

91/15 

198/48 

gg 

26 

150ct71 

1 

5.7 

-14.1 

- 73.3 

268/2 

177/24 

gg 

27 

31 Jan55 

1 

6.75 -12.5 

- 57.4 

142/1 

330/89 

w 

28 

13 Feb 6 4 

1 

5.3 

-18.1 

- 56.8 

211/7 

301/62 

w 

29 

1975-1976 4 

1 

<3.8 

-20.2 

- 44.7 

140/10 


w 

30 


2 


ca 62.8 

ca 6.5 

271 


q 

31 


2 


ca 59.9 

ca 5.5 

86 


q 

32 


2 


ca 5S.1 

ca 7.0 

294 


/“T 
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Table 1 continued 


Number 1 

Date 

Type 2 

Lat°N 

Long°E 

P-Axis/Pl 3 

T-Axis/?1 

Ref. 

33 


2 

ca 52.6 

ca - 8.6 

344 


a 

34 


2 

ca 41. 

ca - 8. 

278 


q 

35 

18Feb71 

1 

5.6 51.0 

5.6 

326/36 

56/3 

a 

36 

24Feb52 

1 

4.9 5 49.3 

8.4 

127/0 

37/46 

c 

37 

21Jun71 

1 

5.7 46.3 

5.7 

315/32 

54/14 

z 

38 

5Feb68 

1 

5.7 46.5 

5.7 

134/7 

314/83 

z 

39 

18Jun68 

1 

4.7 45.7 

8.0 

186/9 

88/16 

b 

40 

15May51 

1 

5.0 45.2 

9.5 

107/33 

202/2 

b 

41 

180ct36 

1 

5.6 46.0 

12.0 

162/24 

64/18 

b 

42 

270ct64 

1 

5.4 47.9 

16.0 

237/4 

329/26 

b 

43 

1946-1965 4 

1 

47.4 

8.6 

119/7 

27/7 

b 

44 

8Feb33 

1 

5.4 48.8 

8.2 

334/7 

72/55 

a 

45 

30Dec35 

1 

5.0 5 48.6 

8.2 

318/20 

70/41 

c 

46 

4Sep59 

1 

4.1 5 48.4 

7.7 

156/42 

57/10 

c 

47 

19Sep65 

1 

4.1 5 48.0 

8.3 

187/32 

86/18 

c 

48 

28Apr61 

1 

4.9 5 47.7 

7.9 

ca 145/10 

ca 45/10 

ee 

49 

25May65 

1 

2.7 5 48.0 

9.8 

0/0 

90/0 

ee 

50 

22Jan70 

1 

4.5 48.4 

9.0 

153/14 

63/0 

b 

51 


2 

ca 45.9 

ca 6.9 

40 


r 

52 


2 

ca 45.9 

ca 8.5 

333 


d 

53 


2 

46.1 

8.8 

140 


m 

54 


2 

46.1 

10.3 

165 


m 

55 


2 

47.1 

8.3 

104 


m 

56 


2 

49.5 

6.4 

126 


m 

57 


2 

49.2 

8.0 

75 


m 

58 


2 

49.0 

8.7 

140 


m 

59 


2 

47.6 

13.1 

160 


m 

60 


2 

49.7 

8.7 

125 


m 

61 


2 

48.2 

9.0 

141 


m 

62 


2 

47.8 

7.6 

165 


m 

63 


4 

ca 48.0 

ca 9. 

300 


dd 

64 


4 

ca 47. 

ca 8. 

125 


1 



114 


Number* 

Date 

Type 2 


Table 1 
Lat°N 

continued 

Long°E 

P-Axis/Pl 3 

T-Axis/Pl 

Ref 

65 

200c t7 2 

1 

5.7 


20.6 

- 29.7 

117/2 

207/2 

aa 

66 

300ct71 

1 

6.0 

- 

0.5 

- 4.8 

147/12 

24/66 

nn 

67 

90ct66 

1 

5.1 


12.6 

30.8 

301/36 

42/15 

nn 

68 

7May64 

1 

6.4 


3.9 

35.1 

65/15 

167/35 

9 

69 

20Mar66 

1 

6.1 


0.6 

30.0 

352/80 

114/6 

11 

70 

23Sep63 

1 

5.8 

- 

16.7 

28.7 

138/75 

283/10 

11 

71 

29Sep69 

1 

5.9 

- 

32.9 

19.7 

90/10 

181/6 

t 

72 


2 



5.8 

- 10.0 

320 


g 

73 


2 


- 

17.3 

14.0 

332 (ave. 

of 4) 

k 

74 


2 


- 

20.3 

30.0 

20 


k 

75 


2 


- 

30.0 

22.2 

342 (ave. 

of 2) 

k 

76 

240ct69 

1 

5.3 


24.8 

72.4 

8/0 

98/55 

e 

77 

10Dec67 

1 

6.0 


17.5 

73.7 

338/36 

65/2 

s 

78 

15Apr64 

1 

5.5 


21.7 

88.0 

181/22 

55/55 

e 

79 

27Mar67 

1 

5.4 


15.6 

80.1 

188/20 

8/70 

e 

80 

100c t70 

1 

5.9 


3.6 

86.2 

172/3 

82/11 

h 

81 

25May64 

1 

5.5 

- 

9.1 

88.9 

311/7 

41/7 

rnm 

82 

26Jun71 

1 

5.8 

- 

5.2 

96.9 

' 320/18 

206/51 

nn 

83 

25Jun74 

1 

6.1 

- 

25.8 

CN 

. 

00 

323/28 

167/60 

hh 

84 

140ct68 

1 

6.0 

- 

31.7 

117.0 

271/13 

20/50 

i 

85 

10Mar70 

1 

5.7 


31.0 

116.5 

282/16 

23/39 

1 

86 

24Mar70 

1 

6.2 

- 

21.9 

126.6 

235/0 

114/61 

i 

87 

2 8 Aug 7 2 

1 

5.6 

- 

25.0 

136.4 

181/1 

237/39 

jj 

88 

27May71 

1 

5.6 

- 

53.9 

150.5 

6 

6 

nn 

89 


2 


- 

20.7 

139.5' 

81 


u 

90 


2 


- 

31.5 

146.0 

97 (ave. 

of 2) 

ii 

91 


2 


- 

42.0 

146. 8 

335 


f 

92 

28Apr68 

1 

5.5 


44.8 

174.5 

70/5 

174/56 

nn 

93 

26Apr73 

1 

5.9 


19.9 

- 155.1 

73/29 

339/8 

PP 

94 

24Sep66 

1 

5.3 


12.0 

- 130.8 

225/1 

161/42 

nn 

95 

3May69 

1 

5.2 


8.3 

- 175.5 

— 6 

6 

nn 

96 

13 Apr 6 7 

1 

5.2 

- 

7.0 

- 151.0 

— 6 

6 

nn 

97 

29 Jul68 

1 

4.9 

- 

7.4 

- 148.2 



6 

nn 
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TABLE 1 

continued 




Number 

l Date 

Type 2 

Lat°N 

Long °E 

P-Axis/P13 

T-Axis/Pl 

Ref, 

98 

27Jun57 

1 

56.3 

117.0 

107/73 

341/17 

y 

99 

5Jan68 

1 

56.5 

121.2 

252/58 

2/12 

y 

100 

14Sep58 

1 

56.7 

121.0 

275/40 

125/6 

y 

101 

29Aug59 

1 

52.7 

107.0 

228/70 

323/2 

y 

102 

18Jan67 

1 

56.7 

120.9 

239/20 

335/15 

y 

103 

14Jun71 

1 

56.2 

123.6 

183/58 

3/32 

y 

104 

4Dec57 

1 

45.3 

99.2 

50/5 

310/55 

y 

105 

23Jun58 

1 

48.7 

102.8 

66/5 

333/34 

y 

106 

5Jan67 

1 

48.2 

102.9 

237/4 

147/0 

y 

107 

20Jan67 

1 

48.1 

103.0 

68/0 

337/64 

y 

108 

4Jul74 

1 

45.1 

94.0 

33/0 

123/0 

y 

109 

29Aug63 

1 

39.6 

74.2 

0/15 

180/75 

y 

110 

HPeb69 

1 

41.4 

79.2 

182/3 

279/57 

y 

111 

14Sep69 

1 

39.6 

74.9 

157/12 

276/78 

y 

112 

5Jun70 

1 

42.5 

78.8 

0/15 

180/75 

y 

113 

29Jul70 

1 

39.9 

77.8 

145/15 

278/62 

y 

114 

23Mar71 

1 

41.5 

79.3 

141/9 

28/72 

y 

115 

15Jun71 

1 

41.5 

79.4 

164/10 

344/80 

y 

116 

29Sep74 

1 

40.4 

78.0 

210/20 

84/60 

y 

117 

26Dec51 

1 

40.0 

95.4 

207/16 

297/4 

oo 

11S 

21May62 

1 

37.1 

95.7 

192/0 

282/65 

oo 

119 

19Apr63 

1 

35.5 

96.4 

54/7 

147/22 

oo 

120 

16Mar64 

1 

37.0 

95.5 

170/33 

354/57 

oo 

121 

24Mar71 

1 

35.5 

98.2 

230/9 

126/38 

oo 

122 

7Mar66 

1 

37.4 

115.0 

69/0 

159/0 

oo 

123 

22Mar66 

1 

37.5 

115.0 

65/0 

155/0 

oo 

124 

22Mar66 

1 

37.7 

115.1 

261/7 

171/7 

oo 

125 

27Mar67 

1 

38.6 

116.6 

67/7 

157/7 

oo 

126 

18Jul69 

1 

38.4 

119.5 

69/12 

335/16 

oo 

127 

5Feb66 

1 

26.2 

103.2 

123/0 

33/7 

oo 

128 

13Feb66 

1 

26.2 

103.2 

119/27 

210/2 

oo 

129 

4jan70 

1 

24.1 

102.5 

346/4 

76/4 

oo 
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Number 

Date 

Type 


Lat°N 

Long°E 

P-Axis/Pl 3 

T-Axis/Pl 

Ref 

130 

HMay64 

r 


28.2 

104.0 

272/3 

182/11 

oo 

131 

9Mar73 

1 

5.5 

-34.0 

150.0 

64/6 

323/65 

X 

132 

3Sep68 

1 

5.5 

20.6 

- 62.2 

299/10 

119/80 

nn 

133 

230ct64 

1 

6.4 

19.8 

- 56.0 

164/3 

258/54 

nn 


* Numbers refer to Figure 7. 

2 

•Type: 1 - fault plane solution; 2 = strain-relief in-situ measurement; 

3 = hydrofracture in-situ measurement? and 4 =* stress-sensitive geologic 
features* 

^ Azimuth and plunge of P-axis for fault plane solutions* Azimuth of 
maximum compressive stress for in-situ measurements and geologic 
features. Azimuth in degrees measured clockwise from north. 

** Composite fault plane solution, 
s 

^ Local magnitude M . 

L 

c 

Fault planes could not be -constrained, but indicate a thrust mechanism. 
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TABLE 3 Shear Stress on the Base 
of the Plates for Model E29 


Absolute Rotation 

Plate Pole 1 Shear Stress 2 , bars 

Lat°N Long°E 


Pacific 

-57.3 

71.4 


2.6 

North 

American 

-44.3 

- 35.5 


2.9 

South 

American 

-58.2 

142.7 


2.7 

Eurasian 

43.2 

- 56.1 

m 

0.6 

African 

66.1 

77.8 


1.4 

Indian 

46.8 

2.0 


5.3 

Antarctic 

11.7 

-148.1 


2.8 

Nazca 

43.0 

- 78.5 


10.1 

Cocos 

38.4 

170.8 


47.0 

Caribbean 

— 

— 


0 

Arabian 

36.1 

46.5 


21.9 

Philippine 

-36.9 

- 24.7 


24.0 


Absolute rotation pole of plate (right-handed rule) with 
respect to presumedly fixed mantle. 

Basal shear stress -D-w x _r at 90° distance from the 
rotation pole, where oj is absolute angular velocity and 
r is radius from the Earth's center. 
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FIGURE CAPTIONS 

Figure 1. Intraplate stress orientation data for the African 
plate. 'Plate boundaries shown as heavy lines. Continental 
boundaries, taken as the thousand fathom bathymetric con- 
tour, shown as light lines. Filled circles, open circles, 
and open triangles, represent stress orientations from 
fault plane solutions, in-situ measurements, and geologic 
structures, respectively. The P-and T- axes for fault plane 
solutions point inward and outward, respectively. Both P- and 
T-axes are shown for strike-slip earthquakes, but only the P-axis 
for predominantly thrust earthquakes and the T-axis for predomi- 
nantly normal faulting earthquakes are indicated. Filled circles 
without arrows denote thrust-fault earthquakes for which 
the orientation of the P-axis is poorly constrained. For in- 
situ data, the line gives the direction of the maximum 
horizontal compressive stress. Dashed lines indicate less 
reliable data. Data for fault plane solutions are from 
Sykes [1967] , Fairhead and Girdler [1971] , Maasha and 
Molnar [1972] , Sykes and Sbar [1974] , and Richardson and 
Solomon [1977] . Data for in-situ stress directions are 
from Gay [1977] and Hast [1969] . 

Figure 2. Intraplate stress orientation data for the Indian 
plate. All conventions as in Figure 1. Data for fault 
plane, solutions are from Sykes [1970] /Fitch [1972] , 

Fitch et al . [1973], Stewart and Denham [1974], Sykes and 
Sbar [1974], Langston [1976], Chandra [1977], Mills and 
Fitch [1977] and Stein and Okal [1978]. Data for in-situ 


stress directions are from Endersbee and Hofto [1963], 
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Mathews and Edwards [1970 ] , and Stephenson and Murray [1970]. 

Figure 3. Intraplate stress orientation data for Europe. All 
conventions as in Figure 1, except that the T-axes for 
those earthquakes with a significant strike-slip component 
are omitted for clarity. Data for fault plane solutions are 
from Schneider et al . [1966], Ahorner et al . [1972] , , Pavoni 
and Peterschmitt [1974], Ahorner and Schneider [1974] and 
Ahorner [1975].. Data for in-situ stress directions are 
from Hast [1958, 1969, 1973] , Capozza et ai . [1971], 

Greiner and lilies [1977] . Data for geological indicators 
of stress are from Greiner [1975] and Scheidegger [1977a, b] . 

Figure 4. Selected intraplate stress orientations inferred from 
earthquake mechanisms in Asia. All conventions as in Figure 1. 
Each data point represents the average orientations of 
P- and T-axes inferred from the mechanisms of at least 
four closely-spaced earthquakes with consistent fault- 
plane solutions [ Molnar et al ., 1973; Tapponnier and 
Molnar , 1977; 1979]. Most of the 

earthquakes indicated are the product of extensive deformation 
of the Eurasian plate along well-developed faults in 
response to the collision of India with Asia [ Molnar et al . , 
1973] , and as such, cannot be strictly interpreted in 
terms of principal intraplate stress directions. 

Figure 5.' Intraplate stress orientations inferred from earth- 
quake mechanisms in the South American plate. All 
conventions as in Figure 1. Because of the density of the 
data, only the P-axis is shown for events near the Nazca 
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plate boundary, although several of the mechanisms have 
a significant strike-slip component. Data from Wagner [1972] , 
Stauder ' [1975] , and Mendiguren and Richter [1978] . 

Figure 6. Selected intraplate stress orientation data for 
North America east of the Intermountain Seismic Belt. 

All conventions as in Figure 1. Data for fault plane 
solutions are from Stauder and Nuttli [1970] , Smith and 
Sbar [1974], Sykes and Sbar [1974], Street et al . [1974], 
Wetmiller [1975] , Hashizume [1977] and Sbar and Sykes 
[1977] . Data for in situ stress directions are from 
Fraser and Pettitt [1962] , Hooker and Johnson [1969] , von 
Schonfeldt [1970], Raleigh et al . [1972], Roegiers [1974], 
Herget [1974], Strubhar et al . [1975], Sbar and Sykes [1977], 
Haimson [1977] and Zoback and Healy [1977] . 

Figure 7. Global summary of intraplate stress orientation data. 

Numbers refer to Table 1. Other conventions as in Figure 1. 
Figure 8. Finite element grid for mid-latitudes, part 1. 

Constant strain triangular elements have dimensions of 5x5x7 
degrees for plate interiors and 3x3x4 degrees for plate 
boundaries. There are a total of 5246 elements for the 
entire grid. 

Figure 9. Finite element grid for mid-latitudes, part 2. 

For other details, see Figure 8. 

Figure 10.. Finite element grid for the north polar region. 

Figure 11. Finite element grid for the south polar region. 

Figure 12. Location of plate boundary elements with lowered 


Young's modulus. 



Figure 13. Schematic illustration of forces potentially 
important for the driving mechanism. 

Figure 14. Forces at ridge boundaries. Forces are applied 
in the direction of the arrows. 

Figure 15. Forces at subduction zone boundaries. Lines 
denote orientation of forces acting toward the trench 
axis from either side of the plate boundary. 

Figure 16. Forces at continental convergence zones. 

Arrows denote orientation of forces. 

Figure 17. Principal horizontal deviatoric stresses in the 
lithosphere for model El ( Table 2). Principal stress 
axes without arrows and with arrows pointing outward 
denote deviator ic compression and tension, respectively. 
Relative magnitude of principal stresses is indicated 
by the length of stress axes. 

Figure. 18. Principal horizontal deviatoric stresses in the 
lithosphere predicted for model E2 (Table 2). For 
other details, see Figure 17. 

Figure 19. Principal horizontal deviatoric stresses in the 
lithosphere for model E4 (Table 2). For other details, 
see Figure 17 . 

Figure 20. Principal horizontal deviatoric stresses in the 
lithosphere for model E17 (Table 2). For other details, 
see Figure 17. 

Figure 21. Principal horizontal deviatoric stresses in the 
lithosphere for model E27 (Table 2). For other details, 
see Figure 17. 



Figure 22. Principal horizontal deviatoric stresses in the 
lithosphere for model E29 ( Thble. 2). For other details, 
see Figure 17. 

Figure 23. Principal horizontal deviatoric stresses in the 
lithosphere for model E31 ( Ihble 2). For other details, 
see Figure 17. 
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III. THREE DIMENSIONAL FINITE ELEMENT MODEL OF TIME DEPENDENT 
DEFORMATION AND STRESS IN THE LITHOSPHERE FOLLOWING EARTHQUAKE 
FAULTING 

Introduction 

We have developed a versatile scheme to model the time 
dependent tectonic deformation and stresses in the Earth. Three 
dimensional finite elements are used to take into account the 
heterogeneities in Earth properties. A unified numerical 
solution approach is incorporated with the wave front 
technique so that an efficient solution can- be implemented 
on an ordinary size computer, and very general rheological 
constitutive laws can be used. Detailed procedures are given in the 
next section. Theoretical model results of the tectonic 
deformation and stress after a thrust fault earthquake and a 
strike slip earthquake are given in the last section. 

Solution Technique 

The central part of the modelling scheme is the finite 
element method. We used the wave front solution approach 
(Irons, 1970) to reduce the problem to a tractable size 
and an efficient solution, and we used a unified numerical 
solution approach (Zienkiewicz and Cormeau, 1974) to deal 
with a wide range of possible rheological constitutive 


laws . 



The wave front approach represents an attempt to 
minimize the in-core storage requirement while giving a 
relatively efficient solution. Gaussian elimination is the 
central process of this technique. However, a global dis- 
placement-force equation is eliminated as soon as it has 
received all of the element stiffness and consistent nodal 
force contributions. The coefficients are moved to outside 
storage and other equations are updated accordingly. So only 
a small portion of stiffness matrix is in core at a time, 
thus achieving a saving in core usage. Assembly and 
elimination processes are not separated. After completion of 
assembly and elimination, the coefficeints are returned to 
core for back substitution in the order opposite to that 
by which they were saved. The wave front solution technique 
for static elastic problems has been implemented at MIT 
and has proved to be quite efficient (Orringer, 1974). 

To deal with different rheological constitutive laws 
and to introduce generality, we divide the total strain e 
into three parts 

£ = e e + E Cp + £ 0 (1) 

0 CP 

where £ is the elastic strain, £° the initial strain and £ ^ 
the creep strain. Quite generally the constitutive law for 
creep strain rate and stress can be put into the form 

£ cp = r a 


( 2 ) 
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where the dot indicates time derivative. T is a symmetric 
matrix (which may depend on stress) and a is the stress. 

The virtual work principle is 

fSc T c d Si - /6u T b d Si - /6u T t d S = 0 (3) 

Si Si Sc" 

where b is the prescirbed body force, t the prescribed boundary 
traction, u the displacement and T indicates transpose. 
Integration is over the volume Si and traction boundary Sc 
respectively. 

Let e=Lu,u=Na (4) 

L is the operator relating strain and displacement, N the 
shape function and a the nodal displacements. Then equation 
(3) becomes 

<$a T [/ (L N) T c d Si - /N T b d Si - /N T t d S] =0 
Si Si" Sc'" 

or /B T c d Si - F = 0 (5) 

where B = L N, F = /N T b d Si + /N T t d S. 

Si" - Sc" 

Using c = D e e = D (e - e c ^ - e ) where D is the set of elastic 
constants, equation (5) can be put into the standard form 


K a = V 


( 6 ) 
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where K = /B T D B d ft 

~ 

and v = F + /B T D e° d ft - /B T D e Cp d ft 

Equation (6) is the set of linear equations to be solved 
using the wave front solution. £ Cp is obtained by stepping 
in time using equation (2) . It can be shown that this 
procedure can be applied to general visco-plastic problems 
(Zienkiewicz and Cormeau, 1974) , including plasticity and 
creep problems as two special cases. Strain hardening and 
the use of non-associate laws do not introduce formulation 
difficulties. Notice also that in later steps, only nodal 
forces are changed. We have only to update the forces in 
the assembly and elimination procedures, and have thus achieved 
an efficient solution. Depending on the configuration of 
the problem, the CPU usage for each later time step is only 
about 5% to 10% of the initial solution. 

Three Dimensional Models for a Thrust Fault Earthquake and 
Strike Slip Earthquake 

Using techniques described in the last section, we cal- 
culated the quasi-static response for two types of earthquakes, 
a thrust fault and a strike-slip fault. In these models, 
we assume that linear viscoelasticity holds. Studies on 
non-linear models are currently underway. 
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Figure 1 shows a schematic diagram of the thrust fault 
model. . The model fault is 400 km in length and dips at 3 0° 
from the surface to a depth of 40 km. The offset between 
opposite sides of the fault is assumed to be 10 meters. The 
assumed rheological properties are given in Table 1. There 
are five elements in each vertical column. The boundary 
conditions are: pinned bottom (700 km depth), free top 

surface, free surfaces at the three sides far from the 
fault, and a reflection boundary condition (zero normal dis- 
placement) for side AA' . 

Figure 2 gives the finite element grid on the free 
surface for half of the model. Due to symmetry, only half 
of the region needs to be modelled. 

Figure 3 shows the horizontal displacement on the 
free surface along the symmetry line (AA' in Figure 2) . 

There is significant time dependent displacement up to one 
thousand km away from the source. 

Figure 4 gives the horizontal plate velocity as a function 
of time at selected - nodes on the symmetry line AA' (indicated 
by circles in Figure 2) . The plate motion is not uniform; 
accelerations continue for decades after the event. 

Figure 5 gives the normal stress perpendicular to fault 
strike at selected locations along the symmetry line AA' 
(indicated by crosses in Figure 2) . The stress decays with 
time near the fault, but increases with time at a distance 
of several hundred kilometers. This is a rather general 



result as long as there is a low viscosity layer beneath the 
lithosphere. Physically, this is easy to understand: the 

dislocation introduces high stress in its vicinity. Due 
to the asthenosphere low viscosity the stresses relax near 
the dislocation while this load diffuses away from the source. 
For this thrust fault earthquake this diffusion effect is not 
very large on the free surface since the fault motion is not 
purely horizontal. We would anticipate a larger diffusion 
effect on the free surface for a strike slip earthquake whose 
fault motion is entirely horizontal. 

The schematic diagram for the strike slip earthquake 
fault is given in Figure 6. The earthquake has a fault length 
of 200 km and extends from the surface to 15 km depth. The 
offset along the fault is assumed to be 10 meters. Again, 
due to symmetry only half of the region needs to be modelled. 
The finite element grid of half of the free surface region is 
given in Figure 7 . The rheological properties are given in 
Table 2. 

Figure 8 gives the horizontal surface displacement 
immediately after the earthquake at nodes for a quadrant 
of the free surface region. Due to symmetry, along the line 
bisecting the fault (AA' in Figure 7) the displacement is 
parallel to the fault; while along the line containing the 
fault (.00' in Figure 7), the displacement is perpendicular 
to the fault except on the fault. Figure 9 and Figure 10 give 
the accumulated horizontal displacement 9.5 and 49 years 
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later, respectively. There is in general an increase in 
displacement with time, along with some change of direction 
(except on symmetry lines) . This increase in displacement 
can be understood as a result of transferring the load from 
the low viscosity layer to the lithosphere as the stresses 
in the low viscosity layer relax. The parallel displacement 
along the symmetry line (OA in Figure 7) is given in 
Figure 11. The detailed time dependence at several locations 
(indicated by circles in Figure 7) is given in Figure 12, 
and velocity as function of time is given in Figure 13. 

Again there is nonuniform plate motion lasting for decades 
after the event. 

Strike slip earthquakes show interesting stress diffusion 
effects in the lithosphere. Figure' 14 shows the time 
dependence of maximum shear stress at several locations near 
the fault (indicated by crosses in Figure 7) . They decrease 
with time monotonically as we expect for high stress areas. 

Away from the fault tip, the stresses drop off with distance 
rapidly, according to St Venant’s Principle, since the 
earthquake dislocation is a balanced force system. However, 
in front of the fault tip, viscoelastic interaction causes 
the stress to increase with time. Figure 15 shows the time 
dependence of maximum shear stress at some locations in front 
of the fault tip (indicated by squares in Figure 7) . At 160 km 
away from the fault center, the stress increases rapidly in 
the first decade. Furthermore, the time dependent stress 



due to earthquakes will reinforce the pre-stress there. The 
obvious conclusion is that the chance of an earthquake 
happening is greater in the area in front of a fault tip. 

The triggered earthquake may again generate its own stress 
perturbation which may cause migration of earthquakes. 
Although this calculation is not meant to model a specific 
region, this migration of earthquakes along a strike-slip 
fault zone is actually observed, for example, along the 
North Anatolian fault (Dewey, 1976; Toksoz et a_l., 1979). 

Summarizing these viscoelastic models: the results 

indicate that great earthquakes cause significant motion at 
distances up to 1,000 km away and that relaxation processes 
last for decades after the event. These conclusions have a 
direct effect on the interpretation of geodetic data. 

The stress diffusion effect can be useful in the analysis of 
earthquake migration in space and time. 
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Table 1 


Thrust Earthquake Model 


Depth 

(km) 

Viscosity 

(poise) 

0 - 

80 

1.0 x 

x,« 

80 - 

180 

1.0 x- 

lo 20 

180 - 

400 

1.0 x 

10 21 

400 - 

700 

1.0 x 

10 2 2 


6 

Young's Modulus = 2.00 x 10° bars, Poisson's ratio = 0.25 


Table 2 

Strike Slip Earthquake Model 


Depth 

(km) 

Viscosity 

(Poise) 

0 - 

30 

1.0 x 

OJ 

o 
1 — 1 

30 - 

180 

1.0 x 

10 20 

180 - 

400 

1.0 x 

10 21 


Young's Modulus = 2.00 x 10 bars, Poisson's ratio = 0.25 
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Fig . 1 . Schematic diagram of thrust fault model 




161 . 



(Mi) dsia 


Fig . 


<s> 

O 

<s> 

<S> 


<s> 


<s> 


CU 


OJ 


i 


3 Horizontal displacement along symmetry line in 
thrust fault model 


-400 





UEL (cm/year) 


163 . 



Fig. 4 Horizontal velocity at several surface nodes along 
symmetry line 
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Fig. 5. Normal stress perpendicular to fault strike 
at locations near symmetry line for thrust 
fault model 


TIME (uear) 



165 . 



slip fault model 








I 200 Km 


167 



6 


s 

o 



OJ 



i 

! c 
J 3 


8 Initial horizontal displacement on free surface 
for strike slip earthquake model 


Fig 


tnodo t SCO . linear scale 



12 00 km 


168 



Fig. 9 Same quantity as figure 8 at 9.5 years later 
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Fig. 13 Horizontal velocity at several locations along the 
symmetry line ^ 
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Fig. 14 Maximum shear stresses near fault 
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